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1 Introduction 

The dilute ultra-cold Bose gas presents a rare opportunity for theoretical physics: it 
has well-characterized interactions, and it is feasible to begin with the full quantum 
theory and subsequently use well-controlled approximations to develop formalisms 
suitable for calculations. These systems can be precisely manipulated and observed 
in experiments and offer a unique chance to compare computational quantum field 
theories directly with experiment. 

Several aspects of experiments present challenges for theory. First, the experiments 
are usually non-equilibrium with long relaxation times and are well-beyond any sort 
of linearized treatment. Second, the harmonic trapping potentials used in experiments 
complicate the traditional many-body methods which are more suitable for uniform 
systems. The low energy collective dynamics and numerous finite- sized aspects of 
this system critically rely on the external potential being treated as a primary consid- 
eration of the theory. 

At zero temperature an almost pure Bose-Einstein condensate (BEC) forms, and 
for a wide range of situations its dynamics are well described by the time-dependent 
Gross-Pitaevskii equation (GPE), e.g., see [53, 189]. This approach assumes that all 
the atoms are well-represented by a single condensate wavefunction, and the GPE 
describes the coherent evolution of this wavefunction neglecting all spontaneous and 
incoherent processes. However, experiments routinely operate in regimes where such 
processes are important and the GPE provides an inadequate physical description, for 
example: 

• At higher temperatures, approaching the condensation temperature, Tc, a sizable 
thermal cloud will be present. Experiments examining collective oscillation fre- 
quencies of BECs found that for temperatures higher than about O.bTc the GPE, 
due to its neglect of the interplay between the condensate and thermal cloud, in- 
correctly predicts the collective mode frequencies and damping [100, 105, 110]. 

• Two nearly-pure BECs colliding produce a halo of atoms scattered onto a spheri- 
cal shell in momentum space [39]. Provided the phase-space density of the scat- 
tered atoms is low, this can be viewed as incoherent scattering of the individual 
atoms in the condensates, and the GPE can be augmented [10] to account for 
these. However, at higher scattered densities Bose stimulation becomes important, 
and a theory which includes both Bose stimulation as well as incoherent scattering 
is required. 

In order to treat these examples and many others it is necessary to formulate a de- 
scription of Bose gases that combines coherent and incoherent physics in a general, 
yet tractable manner. The key to a successful theoretical approach is the recogni- 
tion that in all of these examples, even when there is no BEC, one or many modes 
of the system have an occupation which is much larger than one quantum. The sys- 
tems are then highly Bose-degenerate, and the matter-wave field behaves much Uke a 
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Figure 1 . (color online) Momentum space density (logarithmic) for classical field simulations at various 
temperatures. Emergence of the condensate is visible as a prominent spike at temperatures below Tc- 

classical field. A set of theoretical approaches relying on the existence of significant 
Bose-degeneracy, known generically as c-field methods, provide a comprehensive 
solution to this problem. 

An example simulation demonstrating such a scenario is shown in figure [T] The 
averaged momentum density of a c-field simulation (see section [3]l which describes 
many degenerate modes of a trapped Bose gas is shown for a range of temperatures 
spanning the critical temperature. The condensate is seen to emerge from the broad 
thermal cloud as the temperature decreases below Tc- 

There are two main unifying features of the c-field techniques we present in this 
review. The first is that the modes of the field theoretic description are divided into 
two regions. The precise way this is done depends on the approach, but generically 
we have: 

C: c-field region. This region is of primary importance in the description of the 
system and is so-named because it is simulated using classical stochastic field 
equations. In the theories we develop here this region is chosen to contain not 
only the condensate, but all other highly Bose-degenerate modes. It may also 
contain modes of low occupation in which important dynamics occurs. 

I: Incoherent region. This region consists of the remaining modes which will indi- 
vidually be sparsely occupied (high energy) thermal or vacuum modes. Depend- 
ing on the temperature and the density of states this region may contain a sig- 
nificant or even dominant fraction of the total number of particles in the system. 
However, they have only a weak influence on the dynamics of the c-field region. 
In the techniques we develop in this review the static and dynamical properties 
of this region will be approximated as being incoherent. 

The second common feature to the c-field techniques is that their evolution equa- 
tions are of similar form to the GPE, but with important modifications. This is the 
primary advantage of the formaUsm — its computational tractability, and capability 
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to simulate experimentally realistic parameter regimes. 

This review is organised as follows. We begin in section |2] where we outline the 
background theory relevant to the application of c-field techniques. We then iden- 
tify three separate implementations of the c-field techniques for different physical 
regimes, which are subsequently described in their own sections. 

These techniques ar^ 

(i) Projected Gross-Pitaevskii equation (PGPE) : In all c-field approaches the 
GPE-like evolution is strictly limited to the C region. This is implemented us- 
ing a projection operator, and when this is the sole modification of the GPE we 
refer to the evolution equation as the projected Gross-Pitaevskii equation. 

The PGPE is used to simulate the c-field region as a micro-canonical system, 
i.e. as an isolated system of fixed energy and number, with all couplings to the in- 
coherent region neglected. This approach is valid for high temperatures {T ~ Tc) 
where the energy cutoff is chosen so that all c-field region modes are highly oc- 
cupied, and quantum fluctuations can be neglected. This theory is discussed in 
section |3j along with applications of this formalism to finite temperature phe- 
nomena. 

(ii) Truncated Wigner PGPE (TWPGPE): When there are modes with low occu- 
pation in the c-field region, additional noise terms must be included in the initial 
conditions to model the quantum-mechanical vacuum fluctuations. Inclusion of 
quantum fluctuations cannot be done exactly, but can be well approximated by 
stochastic sampling of a Wigner distribution for the initial state of the system. 
The method introduces spontaneous processes which are absent in the pure GPE 
theory for which all scattering is stimulated. This formalism underlies all the c- 
field techniques and is presented in section [23] with applications of the theory to 
the non-equilibrium dynamics of systems utT <^Tc considered in section |4] 

(iii) Stochastic PGPE (SPGPE): When exchange of energy and matter between the 
c-field region and the incoherent region is important, additional noise terms ap- 
pear in the theory as well as in the initial conditions, via the truncated Wigner 
function method as above. 

This approach is applicable in the same temperature regime as the PGPE how- 
ever it differs from that formalism in that scattering processes, which couple to the 
incoherent region, are included. The theory is implemented by solving the PGPE 
with additional dissipative and stochastic terms. This transforms the description 
of the c-field region to a grand canonical form which includes the exchange of 
particles and energy between the regions. This method, discussed in section [5} is 
well suited for modeling the dynamics of evaporative cooling and, for example, 
vortex formation during Bose-Einstein condensation. 



The term "c-field" techniques has been coined to unify the methods discussed in this review. This terminology 
derives from the name, "classical field method", often given to the pure PGPE formalism, but which we have avoided 
because it can give the misleading impression that is not a quantum mechanical treatment. 
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In all the above c-field techniques it is important to ensure the numerical solutions 
of the equations inside the c-field region do not develop components outside that re- 
gion. Significant research has gone into developing numerical methods for efficiently 
evolving projected equations, particularly the challenge of implementing a projection 
operator efficiently and without compromising the tractability of the equation com- 
pared to the usual GPE. This is discussed further in Appendix [A] and Appendix [B] 
In this review we will show that a wide range of problems can be solved using these 
methods, and that accurate and reliable quantitative results can be computed. 



2 Background formalism 

2.1 Effective field theory for the dilute Bose gas 

In this section we develop the basic formalism for the review. We begin by restricting 
the full Hamiltonian to a low energy subspace, L, for which an effective field theory 
provides an accurate description of the gas with a contact interaction. We then further 
divide this subspace into the C and I regions central to our development of the c-field 
techniques. Our basic approach here follows the derivation given in [74]. 

Our starting point for describing a system of bosonic atoms interacting via an 
interatomic potential U (x) is the second quantized Hamiltonian 

H = J^t/\4''(x)//sp*(x)-F ^ JJ't/^xt/^x't'(x)t^(x0t/(x-x04'(x'mx), (1) 

where ^(x) is the bosonic field operator, and 

H,p^Ho + 6V{x,t), (2) 

Ho = -— + Voix), (3) 
2m 

are the single particle and basis Hamiltonians respectively, with Vq{x) the external 
potential. These Hamiltonians differ by the inclusion of a "perturbation potential" 
6V{x, t), which we include for generality. The basis Hamiltonian, Hq, is so named 
because we will use its eigenstates as a basis for the low-energy description of the 



system, in particular to define the c-field region in section 2.2 The inter-atomic po- 
tential, U{x), has a size characterized by the effective range parameter tq, and only 
depends on the relative separation of the atoms. 

In typical ultra-cold atom experiments the length scales of interest are much greater 
than ro, and the full details of the inter-atomic potential are unnecessary. It is desir- 
able, therefore, to develop a theory that eliminates the need to consider such small 
length scales, and hence the microscopic details of the coUisional interaction can be 
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parameterized in terms of the S-wave scattering length alone - such an approach is 
known as an effective field theory. 

Formally this procedure can be implemented by restricting our attention to a low- 
energy subspace, L, that is spanned by single particle states of energy less than an 
appropriately chosen energy cutoff Smax- This eliminates all momentum states with 
momentum exceeding hA{x) - ^/2m{Em'dx - Vo(x)) at x, and in doing so effectively 
"coarse grains" our description to a length scale of 1/A(x). While our choice of ^max 
is in principle arbitrary, the following criteria ensure a simple and accurate effective 
field theory emerges: 

(i) £max h^/2mr^, so that we eliminate the need to include short wavelength com- 
ponents of the wavefunction that occur in the interaction region. Integrating out 
these high energy states allows the inter-atomic interaction to be replaced by the 
two body J-matrix, which in the zero energy limit becomes [74, 162] 



where Us is the S-wave scattering length, 
(ii) Emax ^ kBT,fi, where // is the chemical potential, so that the eliminated states will 
not be occupied by thermal or interaction effects. This requirement ensures that 
the J-matrix does not depend on the population of states that are eliminated in the 
theory, i.e., avoiding the need to consider a many -body J-matrix. 

As long as these conditions are satisfied, the effective field theory derived should be 
insensitive to the precise value of ^max used. 

We can introduce a coarse-grained field operator, (^(x), which only contains modes 
in L, and is described by the effective Hamiltonian 



It must be emphasized that this resulting field theory has a cutoff, so that the com- 
mutation relations of these new field operators are not precise delta functions: 




X ^' (x)«A^(x)(A(x)«A(x). 



(5) 



f0'(x),^^(x')] =(5l(x-x'). 
In equation ([5]) we have introduced the coupling constant 



(6) 
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(7) 



u = 



(2n)\ 



'm 



where the integral is taken over the momentum space of the L-region and accounts 
for the cutoff dependence of the coupling constant (e.g. see appendix A of [194]). 
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Position 



Figure 2. Schematic view of the c-field region, the incoherent region, and ehminated states for a harmonic trap. 
The c-field atoms require a quantum description, and incoherent atoms may be treated using quantum kinetic theory 



For the special case where the potential is slowly varying compared to the local 
cutoff wavevector A(x), we have 5l(x - x') - sin(A(x)|x' -x|)/27r^|x' -xp. The 
function 5l plays the role of a kind of coarse-grained delta function which in general 
has a spatially dependent width; however, it is also a projector into the subspace of 
non-eliminated modes. Using the commutation relation ([6]) the Heisenberg equation 
of motion for the corresponding field operator takes the form 

ifi^-^ - J d'^ <5l(x - x') {//sp<A(x') + u^\x'mx')!j,{ii!)\ . (8) 

The main purpose of the methods discussed in this review is to simulate this equation 
in various regimes. 



2.2 Projection into the c-field region 

2.2.1 Projection operators 



In section 2.1 we developed an effective field theory description of the cold-atom 
Hamiltonian derived by eliminating states outside of the L region. The resulting ef- 
fective Hamiltonian ([5]) and equation of motion ([8]) are restricted to this space. 
We now turn to a quantitative definition of the L region. This is accomplished by 
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expanding the coarse-grained field operator as 

iff(x) ^ ^ a„</'„(x), (9) 

«eL 

where (pni'^) are single particle eigenstates of the basis Hamiltonian with energy e„, 
i.e. 



eMx) = //o0„(x). (10) 

The operators a„ satisfy the usual Bose commutation relations, [a,, 5/] = 0, and 
[ai,a'j] - 6ij. The restriction of the summation in equation (9 1 to modes in L is 
defined by L = {« : e„ < £max}- 

In general the requirements put on £max for a useful effective field theory to emerge 
lead to an L-space that is far too large to simulate. Furthermore, the validity condi- 
tions for the c-field methods typically restrict their application to describing a sub- 
system of L. Thus it is necessary to further subdivide L into two regions: 

(i) The c-field region (C), which will normally consist of the lowest energy modes 
in L and will be numerically simulated using classical fields{^ 

(ii) The incoherent region (I), consisting of all the modes of L not in C. The choice 
of I will be such that any atoms occupying this region will be best described by a 
particle-like description. 

For all cases considered in this review, these regions are defined by a single particle 
energy eciit> such that C is spanned by the single particle modes with energy e < ^cut 
and I is spanned by the single particle modes with energy ecui < e < £max- We define 
projectors for these regions as 

^PclTO}^ J]<^„(x) r d'x'f,ix')Fix'), (11) 

neC ^ 

Pi{F(x)}^J]Ux) f d\'f„{x')F{x'), (12) 



where C = {n : €„ < fcut) and I = {n : tcut < ^ £^maxl, such that L = C -i- 1 and 
n^c = 0. 



' It is worth noting that the C region has also been refeired to as the coherent region, the condensate band, or the 
classical region in the literature, although some of these names also imply additional restrictions on C. 
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We define quantum field operators for the c-field and incoherent regions as 

<Ac(x) = Pdk^)} - 2 ^«'^«(^^' ^^^^ 

neC 

.Ai(x) = !Pi{(A(x)} = Yj ^"M^l (14) 
)7el 

«A(x) - <Ac(x) + to- (15) 

Most of the theoretical development in this review will be made in terms of the i^c(x) 
and lAi(x) operators. 

Properties of projectors. It is important to note that the kernel of the ;Pc -projector, 
namely 

5c(x,x')^2]<^„(x)C(x'), (16) 

neC 

is the commutator of the c-field operator, i.e. 

[{frc{x),^l{x')]=6c(x,x'). (17) 
The function dc plays the role of a Dirac-delta function for any function in C, e.g.. 



JV<5c(x,x')(Ac(x') = <Ac(x), (18) 



which also follows from the idempotence property, Pc'Pc = "Pc- The imposition 
of an energy cutoff thus has major consequences for the field theory. Imposing a 
cutoff in momentum leads to a discretized form of the continuous field theory with 
6c{Xj,xt) = 6jk/AV for a lattice representation with volume AV per lattice point. 
Imposing it in another basis leads to a field theory with a position dependent com- 
mutator ( 16l. 

2.2.2 The Hamiltonian and equation of motion 

The effective Hamitonian can now be rewritten in terms of i^c and i^i, and because 
the projection is in terms of eigenfunctions of the single-particle Hamiltonian, cross 
terms in the result will appear only in the quartic interaction term. The different ways 
of approaching the implementation of the c-field theory depend on the way in which 
these cross-terms, which connect C and I are dealt with. The three cases are: 

(i) Projected Gross-Pitaevskii equation: The cross terms are dropped, but it is as- 
sumed that all mode occupations in C are significantly greater than one through- 



12 



out the dynamical evolution. Thus, there is significant occupation of I, but this is 
taken as fully thermalized, and conditions are chosen so that when the motion in 
C reaches equilibrium it matches smoothly to I. 

(ii) Truncated Wigner PGPE: Here one deals with processes in which many modes 
in C are unoccupied, including all the higher modes. In this case I is unoccupied, 
and the effect of the cross terms is negligible, and they are dropped. However, the 
quantum fluctuations in C have a significant effect, and this is taken into account 
by including a random element corresponding to half a quantum occupation in 
each mode. 

(iii) Stochastic PGPE : Here one accounts for interactions between C and I by as- 
suming I is thermally occupied, and by using quantum stochastic techniques, 
terms which involve dynamic noise and damping are introduced. When the trun- 
cated Wigner PGPE is also used, the resulting equation of motion is a modified 
PGPE with noise and damping. 

These methods differ in both their conditions for validity and in the details of their 
numerical implementation. We stress that in general a system may evolve from a 
regime where (ii) is the best description (zero temperature BEC), through to a regime 
where (i) is applicable (a thermalised classical field), and finally through sufficient 
heating, to the realm where (iii) may be appropriate. The parameters determining 
this crossover, and therefore the appropriate definition of C and I, are the mode 
occupancies, strength of interactions, and temperature of the system which are in 
general time dependent. Thus, some care must be taken when applying the methods 



and we shall return to this point when discussing validity criteria in section 2.3.9 



2.3 Wigner formalism and the truncated Wigner approximation 



The justification for all three c-field techniques can be made using a Wigner distri- 
bution methodology. Discussion of the properties of the Wigner distribution can be 
found in many places (see for example [78,211]). Here we provide an introduction to 
the theory using the single mode case in sections 2.3.1 and 2.3.2[ before discussing 



the quantum field case in section 2.3.3 



2.3.1 Wigner representation of a single quantum mode 

The Fourier transform of a classical probability distribution is known as its char- 
acteristic function, and moments of the distribution are proportional to derivatives 
of the characteristic function. This same procedure can be generalized to a system 
described by a quantum density operator: For a single bosonic mode with density 
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operator p, the symmetrically orderec^characteristic function is defined as 



Xw(A, A*) = tr[pe''' (19) 

where A is a complex variable. The symmetrically ordered moments are given by the 
derivative ofxw at A = 0, i.e. 

(l^'(^'n,XI;)'(-^)'-<-''^-'l...' 

where {}sym means a symmetrical product of the operators, which is a average of all 
ways of ordering the operators, e.g. 

|aa'| =-|aa'+fl'a|, (21) 
a ^a' j j = - <yya j a + a aa a + a a a + aya j a + aa aa (22) 



The Wigner function was introduced by Wigner in 1932 [213] and is defined as a 
Fourier transform of the symmetrically ordered quantum characteristic function 

W{a, a*) = ^ J (f-A e^'"-^"' xw{A, A*). (23) 

The Wigner function exists for any density matrix (see [78] for a proof), and in 
section!: 



2.3.7 we give the Wigner functions for several standard quantum states. 



Integrating equation ([23]) by parts we see that 



cf-a a\ayW{a,a*) = |-^J xw(A,A%^^. (24) 



Thus (from equation ( 20 1) the moments of the Wigner function give symmetrically 
ordered operator averages 



J (fa a\a*y\ 



{{a\a'Y},y^) = a'{a*y = \ d'-a a\a*YW{a,a*), (25) 



There are several different operator orderings that are commonly used to define the quantum characteristic function, 
with the symmetric case being the standard choice for defining the Wigner function. 
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where we have introduced the notation F(a, a*) for averaging a function of phase 
space variables F(a, a*) over the Wigner distribution. This suggests that the Wigner 
function acts like a probability distribution, indeed W(a, a*) is commonly refereed 
to as a quasi-prohahility since it need not be positive. However, for many impor- 
tant classes of quantum states the Wigner function is either positive (or is well- 
approximated by a positive functi on) and c an be interpreted as a probability distri- 
bution. In these cases the average F{a, a*) is equivalently calculated by statistically 
sampling a as a random variable from this distribution and calculating the average 
of F{a, a*) over many such samples. 

Correlation functions of experimental interest are often normally ordered, requir- 
ing some tedious reordering in order to calculate them from symmetrically ordered 
Wigner averages. However, for a normally ordered operator in the form 

0{a,a^) = Y,Cnma^"a"', (26) 

n,m 

it can be shown that the kernel, Ow{(x, a*), for the equivalent stochastic average over 
the Wigner function, O-wia, a*) = {0(a, a^)) is given by [187] 



o„<«,„-)=2^""<-»"(|4)"(<l*l)" 



z=z'=0 



(27) 



giving, for example 



(a^a) = lap - i, (28) 
<atV) = M^-2|^+^. (29) 



An explicit form for Ow(a, a*) can also be found by evaluating the expression 
Owia, «*) = ^ ^^'7 e~'''''^^0(ar - tj/2, a* + rf 12), 



(30) 



which can be useful in certain circumstances, e.g. for reordering exponential opera- 
tors [167]. For example, choosing 0{a, a') - a^a, we obtain 

Ow{a,a*) ^^jd^^ e-\^\"l\a* + r,y2)(a - r,/2) = |ap - ^. (31) 
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2.3.2 Operator correspondences and equations of motion 

The equation of motion for the density operator under Hamiltonian evolution is von 
Neumann's equation 

= [h,p] , (32) 
where H is the Hamiltonian. For typical Hamiltonians the right hand side of equation 



(32 1 will involve products of operators and the density operator, and here we discuss 
how this equivalently maps onto a differential operator acting on the Wigner function. 
Consider for example the operator product ap. Using 

(e.g. see Baker-Hausdorff formula discussion in [78]) and the invariance of the trace 



under cyclic permutation we have (c.f. equation ( 19 1) 



tr [ape^^'-^""] = tr [pe^"' e-^'^e-^^^''^ (34) 



Fourier transforming equation (35 1 and integrating by parts we obtain the correspon- 
dence 

ap^la+l^]w(a,a*). (36) 



2 da 

Similarly, the other correspondences are 



a+p^(a*-i^Jw(a,a*), (37) 
l_d_ 

1_5 
25a 



pa^[a-~\W{a,a*), (38) 



pa' ^\a* + l-^\Wia,a*). (39) 



We now have a set of mappings from operator equations involving the density matrix 
to a partial differential equation for the Wigner function. 
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Application to the damped and driven harmonic oscillator. As a demonstration of the 
Wigner formalism we consider the driven harmonic oscillator with Hamiltonian 



H = fioja'a + fi{ga' + g*a). 



At 



(40) 



where ca is the oscillator frequency, and g is the driving strength. Such an Hamil- 
tonian arises when considering a single mode of an optical resonator driven by a 
coherent laser field. Damping to a vacuum field outside the cavity adds additional 
non-Hamiltonian terms, leading to the master equation 



dp i ^ y ( J. J. i \ 

— = - - \H, p] + - Lapa -a ap - pa a] . 
ot n 2 ^ ' 



(41) 



Using the operator correspondences (36l-(39l we obtain the equivalent equation of 
motion for the Wigner function 



5W 



y ■\ d ( . ^ y ^ . ^ 
— [lua + 2 ^ '^j 5^ \^iioa + -a - ig 



+ - 



(■■ 

y 52 



W{a, a,t) 



2 dada' 



■W{a,a,t) 



(42) 



This evolution equation is of the form of a Fokker- Planck equation (FPE) with a drift 
(first derivative) term and a diffusion (second derivative) term. Here an important 
tool emerges which is central to the techniques in this review: if the initial Wigner 
distribution W{a, a* , 0) is positive and our interest is in moments of the Wigner dis- 



tribution at some later time (e.g. equation ( 25 1) then we can avoid solving the partial 



differential equation ( 42 1 and instead simulate a large number of trajectories gov- 
erned by the (Ito) stochastic differential equation (SDE) 



da = (-ioj - y/2) adt - igdt + yjy/2dw{t), 



(43) 



with the initial conditions a{0) sampled from W{a,a* ,0). Here the dissipative pro- 
cess has generated a diffusive term in the SDE which is a complex Gaussian noise 
satisfying dw{t) = and dw*{t)dw{t) = dt. The justification for this procedure is the 
standard mapping of a Fokker-Planck equation onto a stochastic differential equa- 
tion [72]. In the case where dissipation is absent, the SDE simplifies to an ordinary 
differential equation. This is the situation for much of the work described in this 
review, but the reasons for the simplification to an ordinary differential equation of 
motion are more subtle and are discussed in detail in section [23^ The SDE ( |43] ) also 
contains the essential technical elements of the phase space mapping required for de- 
riving the SPGPE described in section |5] A more complete discussion of the corre- 
spondence between SDEs and FPEs is presented in Appendix|C] Moving to spatially 
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continuous, modally finite dimensional field theories adds some furtlier complexities 
which are addressed in the following sections. 

2.3.3 Adaption to quantum field theory in the c-field region 

The extension from the single mode case to quantum field theory is accomplished 
using the projected functional generalization of the single mode formalism. Because 
there are only a finite number of modes, the theory can be generalized in a form 
which involves minimal additional calculus. For a system with M modes in the c- 
field region, we define the vector of mode amplitudes a = [ao^ Q-i, . . . , aM-\Y and 
the notation 



d^a„. (44) 



neC 

The multimode Wigner function is then given by 

72 



Wc{a,a*) = j ^exp{A'a-a'A)x^iA,A*), (45) 

where = (a*)^, and Xw is the characteristic function for the c-field region den- 
sity operator, pc- Moments of the Wigner distribution give symmetrically ordered 
operator averages, for example 



d'a laJ'Wcia, a*) = { "^'''^ ^/"""i ) . (46) 



Introducing the c-number c-field (c.f. equation (13 1) 

<AC(X) = Yj (^n'Pni^), (47) 



neC 



the field density average corresponding to equation (46 1 is 



,2 . m2txw *n /'AcW'Ac(x) + (Ac(x)iAc(x)\ 

d^a |.Ac(x)rWc(or,Gr*) = / ^\ , (48) 

/-+ - \ <Jr(x, x) 
= (<Ac(x)«Ac(x)) + -^^^. (49) 



The contribution from the projector in equation ( [49] ) represents a central result of 
the Wigner representation of quantum field theory. Physics beyond mean-field the- 
ory arises in the truncated Wigner method because of vacuum noise evident in the 
commutator term (5c(x, x) ([TT]), which accounts for half a quantum per mode of noise 



18 



present in the theory. This noise mimics the role of vacuum fluctuations, but would 
render the theory ultraviolet divergent if all physically allowed modes were included. 
However, the projection into the c-field region involves a finite number of basis func- 
tions, so that in this situation the term is a well defined finite contribution to the 
stochastic average. 



2.3.4 Functional derivative notation 

There is a useful connection between projectors and functional calculus that greatly 
simplifies multimode calculations while still including all the necessary projectors 
into low energy modes. We define the projected derivative operators as 



S d 



(50) 
(51) 



2. 3. 5 Operator correspondences 

Using the projected functional derivatives, one then finds functional operator corre- 
spondences between the density operator, pc, and the Wigner function [78] 



iAc(x)pc < — > 



'Ac(x) + - 



1 6 



i^^(x)pc < — > Uci^) 



1 6 



2 5(Ac(x) 



16, 

Pc^^c(x)^|^c(x)--^|H^c, 
pc^t(x)^|^*(x)4^|Wc, 



(52) 
(53) 
(54) 
(55) 



which are used to map the equation of motion for the density operator to an equation 
of motion for Wc- 



19 



2.3.6 Truncated Wigner approximation 

From equation ^ we see that the time development of tpc in isolatiorQis governed 
by the Hamiltonian 



^(x)//,p«Ac(x) + ^ J (A^(x)0'^;(x)iAc(x)<Ac(x). (56) 



The equation of motion for the density operator in C is then von Neumann's equation 

(57) 



ih— — = [^c,pc(0]- 



dt 



Using the operator correspondences (52 1 - (55 1, the Hamiltonian (56 1 generates the 
time evolution equation 



dWr 



dt 



/"Mi 



:'Ac(x)^ 



ii^ ^ v4?i5<Ac(x)5(A^(x)"^ 5.Ac(x) 
/ ~6 



+ h.c. 



fi (5(Ac(x) 



(//sp + "[|iAc(x)P - 5c(x,x)])^c(x) + h.c. Wc, 



(58) 



where h.c. represents the Hermitian conjugate. Equation (58 1 as it stands is very dif- 
ficult to solve. However, if we are able to neglect the first line of right hand side 
terms, i.e., those containing third order derivatives, then progress can be made. This 
approximation, which is referred to as the truncated Wigner approximation (TWA), 
is valid over a wide regime for the quantum degenerate gas. The resulting description 
is also obtained formally in the classical limit which we describe below. (We discuss 



the basic validity conditions for this approximation further in section 2.3.9 1. As dis- 
cussed in Appendix |C] a mapping to ordinary stochastic differential equations is not 
possible for equation ( [58] ). However, making the TWA, the Wigner function evolu- 
tion takes the form of a Fokker-Planck equation with drift but no diffusion terms. 



I.e., 



dWr 



dt 



He 



/ ^'^{^^^^(^^P + "[l'Ac(x)l'-(5c(x,x)])«Ac(x) + h.c.|wc. (59) 



The Fokker-Planck evolution can be equivalently mapped to a stochastic partial dif- 
ferential equation [72] that describes the trajectory of a single realisation of the field 
i/rc(x), which we refer to as the truncated Wigner projected Gross-Pitaevskii equation 



'We consider the description of coupled C and I regions in section|5] 



20 



(TWPGPE) 



#c(x) 
dt 



Pc {(^sp + "[l-AcCx)!' - 5c(x,x)])<Ac(x)) . 



(60) 



The lack of a diffusion term in ( 59 1 means that no explicit noise term appears in the 
TWPGPE, however as we shall discuss further in section [2.3.7| the initial conditions 
are stochastic and need to be appropriately sampled from the initial Wigner function. 
We remark that the equation of motion ( [59] ) is also known as a Liouville equation. A 
formal property of the Liouville equation is that the distribution function is constant 
along any trajectory in phase space. This can be seen by applying the method of 
characteristics to ( [59] ), which shows that, within the TWA, the Wigner function is 
constant along classical trajectories given by ([60]). 



Classical limit. While we consider the validity conditions for the truncation in section 
2.3.91 here we show that the truncation is exact in the classical limit, which we define 



as 



0, mA^c - constant. 



(61) 



where A^c is the number of c-field region particles 



^c = J d\ l'Ac(x)|^ (classical limit). (62) 
This expression for A^c is only valid in the classical limit, and the general case. 



obtained from equation (49l, is A^c = J'^^x [|i/'c(x)P - (5c(x,x)/2]. The integral 
Jj^^x (Jc(x, x)/2 = M/2, representing the half quantum per mode vacuum noise 
included in the Wigner description. 



Renormalising the c-field according to 0c = vA^cAc. equation (58 1 becomes 



dt 



He 



+ h.c. 



n <50c(x) 



(H,p + mA^cI0c(x)|>c(x) + h.c.lwc, 



(63) 



In the classical limit we have u/Nq 0, so that the third order derivatives vanish 



in equation (63 1, and the TWPGPE (60) is the asymptotically exact equation of mo- 
tion for the system. However, stochastic initial conditions are still present, reflecting 
quantum or thermal fluctuations, or uncertainties in the initial data of the problem. 

A purely deterministic classical description is recovered when the initial state 
approaches a delta function in phase space which is precisely the limit obtained 
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for a zero temperature BEC in a coherent state. For A^c atoms in a single mode 
coherent state with mode function ^o(x), the renormalised field can be written as 
0c(x) = a^o(x), with phase space distribution 



W(a,a*) = — exp -2A^c 



a - 



2\ 



(64) 



where |aoP = A^c- In the classical limit W{a, a*) 6^^\a - 1), giving the TWPGPE 
for the system dynamics with non-stochastic initial conditions. In general thermal 
fluctuations will be present even in the classical limit, and the problem remains a 
stochastic one, subject to deterministic evolution. We note that in the classical limit 
vacuum noise is unimportant, but it can play an important role in BEC physics where 
A^c ~ 10^-10^. Indeed, the effect of zero point fluctuations can be rather striking and 
even the dominant efi'ect in certain circumstances. An important example is given by 



the dynamics of condensate collisions, described in section 4. 1 



Classical mechanics treatment. For future reference, we note that replacing the field 
operator by the classical field (Ac in Hamiltonian ([56|) yields the Hamiltonian 



Hc = j d'x .Ac(x)//sp«Ac(x) + ^ J d'x |(A^(x)|\ (65) 

which we shall also refer to as the energy functional for the field i/'c- The classical 
equation of motion can then be found by defining the Poisson bracket {F, G) for any 
two functional F and G of the classical field i/'c(x) as 



{F,G} 



6F 6G 



~6F 6G 



(5.Ac(x) (5<Ac(x) <5(Ac(x) 5(Ac(x) 
The equation of motion is then found as 



..#c(x) SHc 
ifi — — = {iAc(x),//c} = 



(66) 



dt 



<5<A^(x) 



(67) 



which yields the projected Gross-Pitaevskii equation 
ih 



"'^"^^ = 'Pc {(H,p + m|(Ac(x)|') ^^c(x)} , 



(68) 



as the classical equation of motion of the system. This equation is equivalent to the 



TWPGPE ( 60 1 in the classical limit where the Sq term can be neglected. The theory 
can be cast in standard Hamiltonian form by introducing the canonical position and 
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momentum variables 



Q„ = \|^J2^lan + al), (69) 
P„ = i^|^{al-an), (70) 



where the a„ are the basis amplitudes of i^c(x) (see (47 1) and e,, are the energies 
of the modes comprising the basis. The Poisson bracket then takes the Hamiltonian 
form and any function F(P„, Qn, t) obeys the equation of motion 



dF \^( dF dHc dF dHA dF 

+ (71) 



dt \dQn dPn dPn dQ„ I dt ' 

We emphasize that the Hamiltonian classical mechanics formulation is not only re- 
covered in the continuum limit, but holds generally as a consequence of including the 
projection operator formally in the definition of the classical field i/'c(x). This Hamil- 



tonian property is used in section 3.2.5 to determine microcanonical thermodynamic 
quantities of the c-field. 

2.3.7 Sampling the Wigner distribution 

We have shown that by making the truncated Wigner approximation, simulations of 
ultra-cold Bose gas dynamics under the Hamiltonian (p6| ) are reduced to sim- 
ulations of the PGPE (or more accurately the TWPGPE ( |60[ )) for an ensemble of 
samples of the initial state of the system. The equation of motion is quite easy to 
solve, but sampling the Wigner distribution for a general many-body system is dif- 
ficult. However, sometimes this sampling issue can be avoided, e.g., in the PGPE 
method a random initial field can be used and allowed to thermalize by evolution 



(see section 3.2.2). 



Bogoliubov formalism. Here our basic aim is to present a procedure for sampling the 
Wigner distribution for a cold (T <s: Tc) Bose condensed cloud in thermal equi- 
librium. In this regime the Bogoliubov method provides an appropriate many-body 
description of the system, provided that the number of non-condensate particles, 
Nex = Nc - No, satisfies A^ex Nc- We briefly review the Bogoliubov formalism, 
and refer to references [38,71,99-101, 144, 176] for a more complete discussion. The 
basic Bogoliubov approach is to expand the field operator in the form 

«Ac(x) = -^^o(x) + y \uiix)bj + v*{x)b]] , (72) 

where is the condensate mode normalized to A^o atoms, {uj{x), Vj{x)] are the quasi- 
particle amplitudes, and {bj, b'j] are quasiparticle operators that satisfy the usual Bose 
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commutation relations. The standard procedure is to take the condensate mode as a 
solution to the time-independent Gross-Pitaevskii equation 

yu^o(x) = //.p^o(x) + u |^o(x)|' ^o(x), (73) 

where // is the chemical potential, and then determine {U j, Vj\ and the quasi-particle 
eigenvalues from the Bogoliubov-de Gennes equations 

ef Uj{x) = [h,p + lu |^o(x)|2 - //] [//X) + M^o(x)' V/x), (74) 
-ef y/x) - [//,p + lu |^o(x)|2 - /i] y/x) + M^*(x)2[//x). (75) 



The expansion in equation ( |72| ) diagonalizes the many-body Hamiltonian ( [56| ) to 
quadratic order in the quasi-particle operators, which is adequate in the regime of 
small condensate depletion, so that the quasiparticle levels are thermally occupied 
according to 

-t-^ 1 

{b'bi) = dn—, , (76) 

= 6ijnj. (77) 

We note that the Bogoliubov modes, {Uj{x), Vj{x)}, are in general not orthogonal to 
the condensate. Even though the correct eigenfrequencies are obtained, orthogonality 
is automatic only for the special case of a uniform system. The correct modes for the 



expansion of the field operator (72i can be recovered from (74i and ( [75] ) by taking 
the projection orthogonal to the condensate [143]. Defining the projector 

n«A(x) - m - A^o '^o(x) J d'x' ^*(x')iA(x') (78) 

the orthogonal modes ai-e given by {m/(x), f,(x)} = {PoUi{x), (P^V-ix))*}. 

Wigner sampling of the Bogoliubov state. By introducing the random variables ao, and 
fi (an M - 1 element vector) in place of the operators uq and {bj}, respectively, the 
Wigner distribution for the Bogoliubov state is appropriately sampled as the stochas- 
tic c-field 

(Ac(x) = ^^o(x) + y [m,(x)/3,- + v*{x)/3*] . (79) 



In the Bogoliubov theory outlined above, the condensate and quasi-particle oc- 
cupations are uncorrelated, i.e. the Wigner distribution is of the separable form 
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(a) (b) (c) 




ai -4 -4 a,- ai -4 -4 -4 -4 



Figure 3. Some possible Wigner functions for tlie condensate mode of a small BEC. (a) Coherent state, (b) 
squeezed state, and (c) number state. Wigner distribution for (Nq) = 10 atoms, and a = a,- + iaj. 



Wc = Woiao, a^)Wqp(fi,fi*), and in the following paragraphs we discuss how these 
can be independently sampled. We note that within a number conserving Bogoli- 
ubov approach additional correlations between the condensate and quasi-particles 
arise [193, 194], providing a better description of the low temperature manybody 
state of the gas. 

Condensate mode: coherent state. To a good approximation the condensate can be re- 
garded as being in a coherent state, for which the Wigner function is 



2 

Wo(a, a*) = - exp (-2\a - aof), 



(80) 



where A^o - lo'oP- For large condensate occupation the finite width of Wo can be 
neglected and for all samples of the condensate amplitude we can take ao * VA^- 

Quasi-particle modes: thermalized states. The quasi-particle levels are thermalized 
modes, with a Wigner distribution of the form of a product of uncorrelated Gaus- 
sian quasi-probability distributions, i.e.. 



Wj(Pj,/3*) - - tanh 



exp 



-2^6/ tanh 



knT 



(81) 
(82) 
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This distribution is sampled by the Gaussian complex random variables, {/3j], with 
properties 

^l=M^i = 0, (83) 

In practice these variables can be generated as 

/3,^ ^«,+ 1/2(^^1 , (85) 

where Xj, yj are normally distributed Gaussian random variates with mean zero and 
unit variance. Sampling the field in this way, we recover the correct symmetrically 
ordered moments, e.g 



- (<Ac(x)>, 
= ^o(x), 
l'Ac(x)P - <{<Ar(x)«Ac(x)}sym>, 



(86) 
(87) 
(88) 



l^o(x)|' + l(<^Ai> + (y6/;))(|M/x)|2 + \vjixf). (89) 



j>0 



Vacuum occupation. We note that even in the zero temperature limit, where rij — > 

0, the random variables /3j have the finite variance lyS^p = 1/2, i.e., each mode of 
the system has on average half an atom of vacuum noise, necessary to ensure the 
symmetrically ordered interpretation of Wigner moments. Thus an attribute of the 
Wigner method is that for a simulation with M modes, M/2 virtual particles (i.e. 
vacuum noise) are included in the field in addition to the A^c real particles. 

2.3.8 Alternative methods for sampling tlie Wigner distribution 

Efficient sampling of a number-conserving Bogoliubov state. Sinatra et al. have shown 
how a number conserving version of the Bogoliubov formalism can be implemented 
via a Brownian motion simulation. This approach, developed in references [191, 193, 
194], has the advantage that it does not require diagonalization of the Bogoliubov de- 
Gennes equations. 



Approximate ground state construction. For a nearly pure condensate the appropriate 
ground state of the system is sampled as the T - limit of expression (79 1. However, 
for many non-equilibrium scenarios, the quasi-particle properties of the low energy 
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modes are unimportant, and the vacuum noise can be added in any basis orthonormal 
to the condensate. That is 

<Ac(x) = ^^o(x) + y l{x)aj, (90) 



where {^/(x)} is an orthonormal basis. The condensate ampUtude, ao, is sampled 



as described below equation (80 1, and the other mode amplitudes are sampled as 



Gaussian random variables with a*^^ = 5,//2. This type of construction is useful in 
collision experiments where the vacuum fluctuation of the high energy modes drive 
scattering events. 

Ideal gas ground state. In the absence of interactions, the ground state Wigner function 
can be sampled as 

^c(x) = 2W«/' (91) 

,/■ 

where 0„(x) are the single particle basis states and the a are sampled according to 



ffo - \Nc and a*aj = (5,//2, for /, j > 0. 



Ideal gas high temperature state. For temperatures above Tc expansion ( 9 1 1 also suf- 
fices to describe the thermahzed state of the system but with all a j sampled as Gaus- 
sian random variables with 

1 



a*aj = 5ij \nBE{ej) + - j , (92) 

where nBEi^j) = {exp[(e,- - fi)/kBT] - 1}"' is the Bose-Einstein distribution. 

More general condensate states. It is possible to consider more general states for the 
condensate, e.g. the number state |A^o)> which has the Wigner function 

Woia, a*) = exp (-2|ff|2)L^„(4|ap), (93) 

where L„(x) is the Laguerre polynomial. The number state Wigner function for A^o = 
10 is shown in figure [3|;. It is non-positive-definite, and is highly oscillatory for 
large numbers which makes exact stochastic sampling difficult. However, for large 
A^o the radial distribution is well approximated by a delta function [73]. A Gaussian 
approximation is thus suitable in this regime and a method for sampling the number 
state Wigner function has been developed and shown to reproduce all moments with 
error ~ 0{1/Nq) [153]. However, more simply we can take ao » VA^^'^> where is 
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a uniformly distributed random phase, e [0, 2n}. Sampling 6 this way preserves the 
f/(l) symmetry of the system. 

Other quantum states, such as squeezed states (see figure [3}5), and crescent states 
can be sampled (e.g. see reference [152]) to investigate their influence on the many- 
body dynamics. 

Adiabatic mapping. In reference [169] Polkovnikov et al. sampled the Wigner func- 
tion of an ideal Bose gas in an optical lattice at T = 0. In the ideal case the bare 
modes, {0y(x)}, form the appropriate basis and the Wigner function can be sampled 



as discussed below equation (91 1. The interactions are slowly ramped up to the de- 
sired value to generate samples of the interacting system (justified by the quantum 
adiabatic theorem). 

2.3.9 Validity criteria for the truncated Wigner method 

The only approximation made in deriving the TWPGPE has been the neglect of third 



order derivatives in the evolution equation for the multimode Wigner function (58 1. 
The complete set of validity conditions for the truncation is still the subject of cur- 
rent research. Polkovnikov, who derived the TWA using a path integral method, 
has obtained expressions for the next order corrections in quantum scattering pro- 
cesses [167, 168] that can be used to assess the validity of any simulation. This ap- 
proach represents a fundamental advance in the formulation and application of the 
truncated Wigner method: it promotes the truncated Wigner method to a controlled 
approximation theory since corrections to the TWA can, at least in principle, be ex- 
plicitly calculated. In practice evaluating the corrections is a challenging task for 
large multi-mode problems, and applications have thus far been restricted to discrete 
lattice systems where the relative strength of interactions to linear system evolution 
is straightforwardly defined. We also make note of comparisons that have been made 
between the TWA and exact results [62, 121, 166] to characterize the limitations of 
the approach for quantum optical systems. 

Several practical conditions for ensuring the reliability of truncated Wigner sim- 
ulations in a variety of regimes have emerged in the literature, and we summarize 
these here. Broadly these conditions fall into two categories: (i) those required to 
ensure consistency of short-time dynamics (relative to the thermalization timescale), 
and (ii) those required for simulations over longer timescales where the system may 
thermaUze. 

Short time evolution: Quantum dynamics. The strict condition of validity of the TW- 
PGPE is that all modes of the c-field region are highly occupied, so that the classical 
limit discussed in section 2.3.6| is approached. In general this condition is rather re- 



strictive, especially for systems well-below the critical temperature. 

Another criterion has been derived by Norrie et al. [148] for factorizable Gaus- 
sian states. In particular, those authors considered the class of states with a Wigner 
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function of the form 



Wc{a, a*) = Y[^-^ exp (-r>, - ajf), (94) 



where the random variable aj, with corresponding orthonormal basis mode ^y(x), 
has mean value a/o, and variance T"'. A sufficient condition for the validity of the 
TWPGPE for these states is 



1 

(iAc(x)iAc(x)> - 2<5c(x,x) 



4 



i.e., the system must have sufficiently high density in position space. Generally, this 
condition is much more readily satisfied than the condition of high mode occupancy. 

Sinatra et al. [194] have made detailed comparisons of the truncated Wigner and 
time-dependent Bogoliubov approaches for a uniform Bose condensate in the regime 
T <s: Tc, where the non-condensate population is much smaller than the condensate. 
They find that the truncated Wigner predictions become inaccurate if the quantum 
noise sampled in the initial condition dominates the number of particles in the sys- 
tem, i.e., the condition for validity is 

M 

y « A^c- (96) 

We note that this condition can be rewritten in terms of the spatial density as 
nix)AV » 1, where n{x) = Nq/V and AV = V/M. Since <5c(x,x) = 1/AV and 



l^jP = l/V for the uniform gas, we see that results (96 1 and (95 1 are equivalent for 
this system. 

Long time evolution: Thermalization. As discussed earlier, when modeling a system 
of M modes an additional half quantum per mode of noise is added on average to 
the initial condition. This introduces M/2 virtual particles into the TWPGPE simu- 
lation which should be subtracted to recover the correct operator averages. However, 
under evolution these virtual particles may thermalize, and change the equilibrium 
properties of the system. 

Sinatra et al. [194] have proposed the condition 

\T - Teiassl « T, (97) 



to ensure the long-time validity of a truncated Wigner calculation for a system, where 
T is the initial temperature of the system, and Tdass is the temperature once the noise 
has thermalized. In practice they find that equation ([97]) is best ensured by limiting 
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Figure 4. Schematic view of the modes in a critical system showing the appropriate choice of fcut- 

the number of modes in the numerical calculation, and that an acceptable description 
is obtained if the largest single particle (quasi-particle) energy in the calculation is 
no more than a few ksT . 

Interactions and dimensionality. It is clear that the conditions listed above cannot be 
considered complete as they do not expUcitly involve interactions or dimensionality. 
In regimes where interactions dominate the system can evolve into a strongly cor- 
related state. Important examples for ultra-cold atoms are the Mott-insulating and 
Bose-glass states that emerge when a system of repulsive bosonic atoms is confined 
in a deep optical lattice [35, 108]. Generally speaking the truncated Wigner approxi- 
mation does not provide a good representation for such strongly correlated systems. 

Another regime in which interactions play an important role is in the critical region 
where the system undergoes strong fluctuations. These critical fluctuations are clas- 
sical in nature and thus amenable to the truncated Wigner treatment, however care is 
needed in choosing C to describe this regime. Typically strong fluctuations occur in 
the infra-red modes up to the energy scale un where n is the density, above which the 
modes are well-described by mean-field theory (e.g. see the discussion in [1 19, 172]). 
To provide an accurate description the c-field region must be chosen to include these 
modes, i.e. we must have ^cut > un, as shown schematically in figure [4] 

The occurrence of critical fluctuations can be identified by the Ginzburg cri- 
teria (e.g. see [67]), which predicts that fluctuations are important for the three- 
dimensional Bose gas only in a narrow temperature range about Tc- Outside of this 
range a pure mean-field description provides a good description of the equilibrium 
system. In contrast, for low dimensional systems strong fluctuations prevail over a 
broad temperature range, and inhibit the formation of phase coherence according to 
the Mermin-Wagner-Hohenberg theorem [94, 140]. For these systems a pure mean- 
field analysis is of limited use, yet a c-field description of the low energy modes (C 
region) with a mean-field description of the high energy states (I region) can be used 



to provide a comprehensive description (see [172-175] and section 3.4.2 1. 

Low dimensional systems exhibit several additional properties that make them 
well-suited to simulation by c-field techniques. First, in lower dimensionality the 



30 



rate of thermalization is significantly reduced witli respect to the tliree dimensional 
case. This suggests that the thermalization of the vacuum noise that occurs in trun- 
cated Wigner evolution will happen more slowly and simulations in the low temper- 
ature regime should be valid over longer time scales. Second, the density of states 
increases more slowly with energy as the dimensionality of the system decreases, 
and thus for ID and 2D systems at finite temperature a larger portion of atoms reside 
in low energy, highly occupied modes. 

In general the study of ID and 2D Bose gases is an active area of research, and 
many aspects of their behaviour, particularly dynamics, are not well-understood. 
While the c-field techniques are widely applicable to describing these systems, a 
significant challenge remains to develop techniques for sampling the Wigner distri- 
bution when the Bogoliubov description fails. A few procedures have been developed 
for quasi- ID lattice systems at T « in [169], and for finite temperature quasi-2D 
trapped gases in [190]. 



Role of projection in ensuring validity. Both the short-time and long-time validity con- 
ditions are sensitive to the number of modes in the c-field region and the maximum 
energy of modes represented. For this reason it is essential to have a projector in the 
formal theory to exert as much control as possible over the modes retained in the 
c-field description. 

A natural question that arises is how dependent are the results of a simulation on 
the energy cutoff we use to define the C and I regions? For the case of a system in the 
critical regime it is clear that a lower bound for ecut is provided by the energy scale un 
as discussed above (see figure|4]). More generally, this lower bound is also appropriate 
for inhomogeneous systems outside of the critical regime. This is because we have 
chosen to implement the projector in the single particle basis (see section 2.2 1, which 
provides a good representation of the low energy modes and a well-defined energy 
cutoff only when ecut > un. 

The criteria for setting an upper bound to ecut varies greatly according to applica- 
tion. A large value of ^cut can be used for short time simulations, subject to the con- 
dition given in equation ( [96| ). For equilibrium situations (or long time simulations 
where relaxation occurs) typically the condition ec^ ~ kgT [194] ensures that all the 
modes in the c-field region are appreciably occupied, and thus accurately described 
by the Wigner approach (we discuss this case in more detail in section [3]l. 

Within the aforementioned guides for choosing the cutoff, there is still an apprecia- 
ble degree of freedom. For example, consider the situation shown in figure |4j where 
a particular choice of tcut is indicated that splits our system of interest into the C 
and I regions, which we take to be described by c-field and mean-field descriptions 
respectively. In this case it is clear that an appreciable number of modes appropri- 
ately described by mean-field theory are included in the C region. Moderate shifts in 
the value of ^cut result in the transfer of some modes, well-described by mean-field 
theory, between the regions. Clearly this scenario will have no effect on the physical 
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predictions for the complete system as the modes shifted between the regions are 
equally well described by both formalisms. However, large changes in tcut will lead 
to problems: If fcut is set too low, then strongly fluctuating modes will inappropri- 
ately appear in the I region. If e^cut is too high then many sparsely occupied modes 
will be included in C. 

While the theoretical motivation for choosing the cutoff is clear, there are only a 
few studies that have examined the dependence of simulation results on the cutoff. 
Of most note: 

• Sinatra et al. [194] studied the damping rates for coherent excitations as the num- 
ber of C modes was varied by an order of magnitude. They found that the damping 
rate varied by a factor of two over this range, and the best results (i.e. those agree- 
ing with the Beliaev-Landau damping result) were obtained for low cutoffs. 

• Bradley et al. [30] derived Ehrenfest relations for a c-field system, which display 
an explicit dependence on the projector. They also presented results showing that 
the macroscopic properties of density distribution and the condensate fraction var- 
ied appreciably when the numerical method used was changed from a spectral ap- 
proach using oscillator states (i.e. implementing the C region using a well-defined 
energy cutoff in the single particle basis) to a uniform grid (i.e. implicit projection 
that only provides a momentum cutoff at the inverse of the grid point spacing). 

2.3.10 Features and interpretation of the truncated Wigner method 

Single trajectory interpretation. Phase space methods provide a practical means for 
calculating correlation functions, which can only be compared with the equivalent 
quantities calculated as an ensemble average of experimental measurements. How- 
ever, for highly occupied fields, the behaviour observed in each trajectory of the 
TWA seems to be typical of that seen in single realizations of experiments, (e.g. see 



the results in section 4.1 and section [43] ). Thus it is plausible that single realizations 
of Wigner trajectories should approximately correspond to a possible outcome of a 
given experiment. This is no more surprising than the observation that the GPE is a 
remarkably good predictor of the dynamics of individual experimental runs, and fol- 
lows from taking the classical limit of the full quantum evolution, with the important 
addition of initial fluctuations. 



Restriction to states with positive Wigner function. However, this is not the full story as 
the Wigner function can be negative for some states (see figure[3|:). Negative Wigner 
functions are not amenable to exact treatment by diffusive processes and so there are 
in fact certain quantum states that are inaccessible to the stochastic sampling methods 
described in this review. There are actually two points here. Firstly, negative Wigner 
functions aie difficult to sample as an initial condition. Secondly, a positive Wigner 
function will not become negative under diffusive TWA evolution. While it may pre- 
dict the mean fields accurately, it may not (and cannot) give the conect correlation 
functions for some processes. These are however, fairly rare with EEC. This restric- 
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tion limits the range of quantum phenomena to states with positive Wigner functions, 
ruling out superpositions of number states and demonstrations of the non-locality of 
quantum mechanics (violation of a Bell inequality). 

Spontaneous scattering. The GPE is fundamentally a theory of stimulated (Bose- 
enhanced) scattering, which does not include spontaneous processes. In particular, 
scattering into initially unoccupied modes will not occur, although this may eventu- 
ally occur in computational simulations due to the gradual accumulation of numerical 
errors. The Wigner method, however, sets an irreducible level of initial fluctuations 
in all modes of the c-field i.e. half an atom of vacuum fluctuations. In efi'ect, sponta- 
neous scattering becomes modeled by weakly seeded stimulated scattering. 

MuUimode averaging. The c-field used in truncated Wigner simulations usually con- 
sists of some large number of stochastically sampled modes M. Many observables of 
interest (e.g. column density, cloud rms-width) depend on the values of a significant 
portion of these modes, so that the statistical fluctuations in the value of such ob- 
servables can be quite small. Often the behaviour and evolution of these observables 
exhibit little difl'erence between independent trajectories. 

Long time dynamics. Sampling the initial state introduces fictitious population into the 
system, i.e. the vacuum noise. In ensemble averaged calculations this is subtracted 
when constructing operator averages from trajectory averages. In single trajectory 
dynamics of the truncated Wigner method, the extra population becomes dynam- 
ically thermalized and indistinguishable from the rest of the field. There are two 
eff'ects at work here. First, the truncation of the equation of motion means that quan- 
tum mechanical corrections, which prevent this thermalization, have been neglected. 
Second, by considering single trajectories, the formal correspondence to operator av- 
erages is lost and the results must be interpreted within the context of classical field 
theory. For long times, the advantage of the truncated Wigner method is that it pro- 
vides a more complete physical picture of the system evolution than the GPE, but it 
must be interpreted with some care. The primary gain is the inclusion of spontaneous 
efl'ects in the dynamics from the outset. 



3 The projected Gross-Pitaevskii equation 

The projected Gross-Pitaevskii equation (PGPE) formalism is valid for degenerate 
Bose gases at finite temperature, a system for which many excited modes (in addition 
to the condensate) of the atomic field may have a high mean occupation, i.e., satis- 
fying the criterion (a^a,) » 1. In this section we describe the PGPE formalism, and 
show how it can be formulated to make quantitative comparisons with experiments 
in this regime. Unlike the other c-field techniques described in this review, the PGPE 
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formalism can be described as a "classical field theory" in the sense of the classical 
limit taken in section l2!3.6[ 



3.1 Classical field description of thermal Bose fields 

The suggestion that the Gross-Pitaevskii equation could be used to describe the dy- 
namical evolution of the Bose field in the limit of large mode occupation was first 
made by Svistunov in 1991 [205], and later by Kagan et al. in references [116-118]. 
Damle et al. [44] first investigated this using numerical calculations in 1996. They 
used the homogeneous GPE with a very weak nonlinearity to study the phase- 
ordering kinetics of a Bose gas in two and three dimensions on small grids, and 
performed a scaling analysis of the growth of the condensate fraction in a tempera- 
ture quench. 

Subsequently Marshall et al. [138] studied equilibration of a harmonically trapped 
Bose gas in 2D using the GPE. They observed changes in the population distribution 
of the bare harmonic oscillator states, and relaxation of the density profile from an 
initial asymmetric form to a radially symmetric one, and interpreted these changes 
as thermalization. 

The introduction of a projection operator to restrict the modes represented by the 
GPE was first reported by Davis et al. [51] for the case of a 3D homogenous gas. At 
a similar time thermalization for a homogeneous system was demonstrated by Goral 
et al. [85], who solved the equations of motion for the mode amplitudes explicitly in 
the classical approximation. 

3.1.1 Importance of the projector and numerical methods 

It has been long known that applying classical field theory to the electromagnetic 
field results in the ultra-violet (UV) catastrophe in which an infinite number of modes 
each have the equipartition share of energy, ksT . Thus it would seem that the effects 
of a UV catastrophe would also impact a classical field description of the Bose gas. 
However, the manifestation of the catastrophe is rather different. The GPE is the 
equation of motion for a classical microcanonical field in which the total energy and 
particle number (field normalization) is conserved. In thermal equilibrium this en- 
ergy is shared equally (equipartitioned) between all system modes. Any numerical 
solution of the GPE is constructed from some finite basis, e.g. choosing an equally 
spaced grid (equivalent to choosing a basis of planewaves in the first Brillouin zone). 
Increasing the number of grid points on which the thermally equilibrated GPE solu- 
tion is constructed means the fixed energy is now shared between a larger number of 
modes, so that the average energy per mode (i.e. temperature of the system) will de- 
crease. In this sense, the results of calculations are dependent on the numerical basis, 
or more correctly the number and nature of the modes contained in the calculation. 



However, as we discuss in section 3.1.2[ if the modes of the calculation correspond to 



only the highly occupied modes of the physical system, then a quantitative descrip- 
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tion of the system can be made. For this reason the use of an expUcit projector in the 



PGPE equation (68 1 is of great importance, because it precisely defines the calcula- 
tion without reference to the numerical implementation. In addition, great care must 
be taken in implementing numerical methods for propagating the PGPE so that all 
modes of the c-field region are evolved accurately in order to avoid spurious dynam- 
ics which can lead to an incorrect representation of the physical system of interest. 
We note that a classical field formalism has been developed using unprojected grid 
methods, summarized in Brewczyk et al. [34]. While such an approach seems suit- 
able for investigating qualitative behaviour of Bose gases in various regimes, it has 
not been applied to quantitative comparison with experiments. 

Use of grid methods. Grid methods are ubiquitous in the solution of the GPE, but care 
must be taken in using these methods as the basis of classical field simulations. For 
example, the cubic nonlinearity in the PGPE can generate momentum components 
up to three times larger than those present in the classical field. In a grid representa- 
tion of the field this leads to aliasing, which corresponds to (unphysical) collisions 
between modes that do not conserve momentum. Grid methods can be used effec- 
tively for numerically solving the PGPE for a uniform gas, if several adjustments 
are made: (i) The projector needs to be explicitly implemented. E.g. the single par- 
ticle energy cutoff, discussed in section |2.2[ can be implemented by setting to zero 
all modes outside a sphere of radius fikcui = yfhne^i in momentum space, (ii) A 
large number of momentum states outside the projected region need to be retained 
to avoid the aliasing problem. The equivalent position space requirement is that if 
^cut is the largest wavevector retained by the projector, then a spatial grid of spacing 
A;c = nllkcui needs to be used to evaluate the nonlinear term, which is twice as dense 
as the Nyquist sampling requirement, ^^xn = n/kcut- Additional discussion of these 
issues is given in Appendix [B] 

The experimentally relevant harmonically trapped system poses a more formidable 
challenge since the natural modes of grid representation (i.e. planewave modes) 
bear little resemblance to the harmonic oscillator modes, making projection difficult. 
Also, we note that in typical experimental regimes there should be of order 10^ - 10"^ 



modes in the c-field region (see section 3.1.2 1, whereas grid methods usually require 
> 10^ points to accurately simulate the GPE in 3D. Details of an efficient numerical 
algorithm for the PGPE in a harmonic trap is summarized in Appendix [A] 



3.1.2 c-field region for tlie PGPE: the "classical region" 

For the PGPE formalism, the occupations of all the modes of the c-field region satisfy 
the mean occupation requirement 



<aja„) w a*an > ricm, (98) 
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where Wcut is a number of order one (typical choices range from about 1 to 10). 
Loosely, ?iciit is the degree of coherence of the mode, and should be compared 
to the basic level of quantum fluctuations set by the Wigner function requirement 
{a*a„)„ac - 5 for vacuum modes. We hence refer to the modes satisfying this con- 
dition as constituting the classical region, in the sense that quantum corrections to 



c-field equation of motion (68 1 for these modes are small (see section 2.3.6 1. A fur- 



10' 




■B ^o' 

CO 
Q. 

O 
O 

o 
c 

CO 

10 



10" 




10" 



Bose-Einstein 



10"' 

(E-tl)/kgT 



10" 



Figure 5. (color online). Comparison of the quantum Bose-Einstein and classical equipartition predictions for the 

mean occupation of a single particle mode. 



ther justification for setting «cut can be obtained from figure [5} where the mean mode 
occupation is examined as a function of the scaled single particle energy. When the 
parameter (e - ii)/kBT <sc 1, then the exponential in the quantum Bose-Einstein dis- 
tribution 



1 



explie - iu)/kBT] - l' 



(99) 



can be expanded to first order to give the classical equipartition distribution 

IcbT 



HEQ{e) 



(100) 



In figure[5]it can be seen that these two distributions are in good agreement for highly 
occupied modes, i.e. modes satisfying e - fi < kgT with mean occupation n> \. 

The properties and size of the classical region for typical experimental parameters 
are not a priori obvious, especially for the case of an interacting gas. In figure |6] 
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we consider the case of a harmonically trapped system and estimate the number of 
classical region modes and the number of particles occupying those modes using a 
Hartree-Fock mean-field calculation [28]. Those results reveal that the number of 
classical modes is maximum at the condensation transition, with of order several 
thousand modes satisfying condition ([98]) for the parameters of this calculation. 



Strictly, equation (98 1 is only applicable to the noninteracting modes of the gas, but 
can be generalized to the interacting system, e.g. by analysing the occupation of the 
natural orbitals of the one-body density matrix. However, typically 6cut is sufficiently 
large compared to the interaction energy scale, that the highest energy noninteracting 
modes (i.e. (p„{x) with e„ w fcut) in C are a good approximation to the modes of the 



interacting system. In this case equation (98 1 can be directly applied to these high 
energy modes. 

3.1.3 PGPE fo rma lism 



The appropriate equation of motion for the c-field is the PGPE (68 1. Since all the 
modes in C are highly occupied, in a three-dimensional system we find an appre- 
ciable number of atoms residing in the incoherent region. Thus, the detailed non- 
equilibrium dynamics of the system will in general depend on the exchange of energy 
and particles between C and I. A consistent formalism for including these processes 
is described in section [5] However, for the purposes of this section, we will assume 
that for near equilibrium scenarios the C and I regions are weakly coupled, and we 
can treat each region in isolation as long as we match their temperatures and chemi- 
cal potentials. Building on this assumption, in the remainder of sectionj3]we mainly 
concern ourselves with the properties and interpretation of the PGPE ([68f as a micro- 
canonical means for describing the classical region of a finite temperature Bose gas. 



In section 3.2.6 we discuss a mean-field treatment of the incoherent region, I, as a 



means to providing a quantitative description of the full system. 



3.2 Hands-on introduction to the PGPE formalism 

In this section we introduce the basic ideas of the PGPE formalism by example. To do 
this we present a hands-on case study of how to prepare initial states and evolve them 
to thermal equilibrium. We introduce various tools for analyzing PGPE simulations 
as needed. 

3.2.1 Simulation parameters 

To guide our presentation we illustrate the PGPE method using simulations for an 
experimentally realistic system. We take this system to be a harmonically trapped 
gas of rubidium-87 atoms in a potential 



(101) 
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Figure 6. (color online). Size and population of the classical region for a harmonically trapped system, (a) 
Classical region population including (line) and excluding (dashed) the condensate occupation, (b) Number of 
classical region modes. Results: pink (grey) N = 2x 10^ atoms, and blue (black) N = I X 10^ atoms. Results 
calculated using Hartree Fock theory (see [28]) for rubidium-87 atoms in an isotropic harmonic trap of frequency 
100 Hz with Hcut = 3 used to define the classical region 



where {cojt, lOy, oj^} are the oscillation frequencies (and with dV = 0) . For definiteness 
we take a)x = 2nx 120 Hz, ojy^^ =27rx30 Hz, i.e. the trap has afat pancake geometry 
with the A:-direction being tightly confined. For the calculations we fix the cutoff 
defining the c-field region at ecut = 337j(jjj., so that there are M = 1560 single particle 
modes in C. We take the number of atoms in this region to be fixed at A^c = 10^ (see 
equation (|62])), and will verify a posteriori that all the modes are highly occupied as 
required for the validity of the PGPE formalism. 

3.2.2 Initial state preparation 

The nonlinearity of the PGPE causes its evolution to be ergodic and many of the is- 
sues involved with appropriate sampling of initial conditions in the truncated Wigner 
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approach can be avoided if we are only interested in the equilibrium properties of the 
system. 

Thus the generic method to study finite temperature regimes is to begin with a 
randomized initial state with some definite energy, as specified by the c-field Hamil- 



tonian Eq = Hc[ipc] <65 1. The c-field energy is a constant of motion for the PGPE 
(68 1 and forms a convenient macroscopic constraint for specifying the thermal state 
of the system. The procedure for making such energy states is rather arbitrary. We 
choose to make use of the Thomas-Fermi approximation to the condensate mode 



^Tf(x) = / ^TF Voix) ^ ^^^^ _ ^^^^^^^ ^ ^^^2) 



where 9{x) is the unit step function and 



is the Thomas-Fermi chemical potential [43], with with 0? = OL)x(^y(^z and a - 
ylfilmoj. We can generate a state of desired energy by superimposing the Thomas- 
Fermi state with a (high energy) randomized state, ^, (x), according to 

^£(x) = Mtf(x) + Pi^.(x), (104) 

where the variables pQ and p\ are adjusted to obtain the desired energy. In practice, 
is approximately orthogonal to ^tf and we can take pi = yjl - \po\^- 
For reference, under the constraint of fixed c-field normalization, A^c> the minimum 
energy configuration corresponds to the zero temperature case with all atoms residing 
in the condensate mode ^o(x) (i.e. A^o - ^c)> which can be obtained by solving the 



To 

we 



time-independent Gross-Pitaevskii equation (73 1. 
3.2.3 PGPE thermalization 

Here we present evidence for thermalization of the PGPE in the trap|)ed system, 
do this we consider two initial microstates of the c-field: ^^'(x) and ^^(x), which 
will refer to as cases (a) and (b) respectively. Both of these initial states have the same 
energy Eq = lO.ONcho)^, and are constructed according to the procedure outlined 
in section 3.2.2 but with a modified choice of the Thomas-Fermi state, ^tf» as we 
discuss below. Such initial states will not be equilibrium states, and during PGPE 
evolution will thermalize. To emphasize the initial non-equilibrium dynamics and 
the role of thermalization we choose to use distorted Thomas-Fermi states in ( |104| ): 



initial state ^^'(x) is produced using a Thomas-Fermi state that has been squeezed 
in the x direction; initial state ^^'(x) is produced using a Thomas-Fermi state that 
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Figure 7. (color online) Relaxation to equilibrium. Density slices of the two non-equilibrium initial states: (x) 

(al)-(a2), and '(x) (bl)-(b2), and the respective states they evolve to at / = lOOOmi (a3)-(a4) and (b3)-(b4). Both 
states have Ec = lO.ONcfia)-. The rms-width of the c-field in the jr-direction (c) during the first few trap periods 

and (d) after 15 trap periods. Width of simulation "case(a) " (black line) and the "case (b)" (grey line). Note: other 
parameters given in section [3.2.1| and colourmap in density plots corresponds to logjQ of the density measured in 

units of ijimy^ . 



has been squeezed in the y direction. These two initial states, while having the same 
energy, are clearly very distinct in spatial character as revealed by the density slices 
shown in figure [7Jal)-(a2) and (bl)-(b2). The final states after PGPE evolution for 
1000 ms are shown in figure [7Ja3)-(a4) and (b3)-(b4). 

These results show that the system thermalizes, in the sense that the system evolves 
to more-likely microstates. Indeed, while the states in figures |7ja3)-(a4) and figures 
|7jb3)-(b4) are not identical (differ by fluctuations), they are much more similar than 
their initial states. 

To examine the dynamics of thermalization more carefully we show the .x-widths 
of the states, as characterized by the rms value x^s = {x^)t - {x}j (where ( }t is the 
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expectation at time t), in figure [Vjc) and (d). Initially, systems (a) and (b) exhibit 
large oscillations and differ strongly in their width dynamics (see Fig|7jc)), reflect- 
ing the differences in the intial non-equilibrium states. After approximately 20 ms 
the large scale width oscillations have damped significantly, leaving much smaller 
fluctuations. In figure[7Jd) we show the width dynamics from 50 ms to 500 ms. Here 
the width fluctuations are about an order of magnitude smaller than the initial os- 
cillations, and despite both systems beginning from very distinct initial states, these 
dynamics are consistent with both systems thermalizing to the same equilibrium, 
i.e. the same mean width and fluctuation properties. 

There are a large variety of observables that we could compute to examine the 
thermalization of the system. However in general we typically find that the system 
relaxes towards equilibrium appreciably within a few trap periods. For typical simu- 
lations, where we are interested in equilibrium properties and start from the (undis- 
torted) initial state described in section |3.2.2[ we evolve the c-field for several tens 
of trap periods to thermalize before sampling for system properties. 

3.2.4 Equilibrium: Ergodicity, correlation functions and condensate fraction 
Ergodic averaging correlation functions. The c-field energy, given by the functional 



(65 1, is a constant of motion under PGPE evolution. Indeed, the energy and other 
such constants of motion, e.g. field normalization and angular momenta (when the 
trap has rotational symmetry), take the form of the macroscopic constraints on the 
thermal state of the system. In principle, the equilibrium properties of the system 
could be determined by ensemble averaging over all fields consistent with these con- 



straints. The nonlinearity of equation (65 1 makes finding all functions ipc, for a given 
normalization and energy, impossible without approximation. If we were to move 
beyond the microcanonical ensemble, some form of Monte Carlo sampling could be 
used, although we do not pursue this possibility here. 

In contrast, numerical methods for evolving the PGPE are well-developed and 
allow a different means to sample the ensemble: we can make use of the ergodic 
hypothesis (that a system will in time visit every accessible configuration in phase 
space without bias) to sample microstates of the system. Thus, for an ergodic system, 
an ensemble average of an observable, O, can be calculated by a time average over a 
sufficiently long period of dynamical evolution, i.e. 

where {f^} e {9i,9 + 0,} is a set of Ms time instances at which the system evolution 
has been sampled [52, 86]. For this choice to be a accurate estimates of the ensemble 
average we require Ms » 1 , and the time span over which averaging is done to be 
long compared to the slowest time scale in the problem, e.g. the longest harmonic 
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oscillator period. 

In general the observables of interest are of the form of a correlation function of 
the field, typically a product of quantum field operators such as 



{6) = (tAc(xi) . . . iAc(x./)iAc(Xj+i) . . . iAc(x«)>. 



(106) 



This expression could also be generalized to multi-time correlations, although we 
will not do so here. Note that here the correlation functions only involve the c-field 
opera tor: We discuss corre lations involv ing i ncoherent region operators in section 
3.2.6 and in section 3.4.3 To evaluate (106 1 we make the substitution i^cCx) 



i/'c(x), transforming the expression to the general classical field form, and then re- 



place the ensemble average with a time-average according to ( 105 1. We note that this 
procedure is in accordance with that outlined for the truncated Wigner approach (see 
2.3. 1| ) as we can neglect the commutation relations of the operator fields for 



section 



the highly occupied modes described by the PGPE. 



e 




Figure 8. (color online) c-field position (a)-(c) and momentum (d)-(f) density. Simulation parameters: (a), (d) 
£c = 24A?cS6j;, (b), (e) Eq = llNcno)-^, (c),(t),(g) Ec = ISNcTioj-. Other parameters: A'c = 10** **^Rb atoms, 
simulations are evolved for / = lO^n/aj,, with Ms = 2500 samples taken over the last half of the simulation used to 
time-average. Trap parameters and cutoff are given in section [3.2.1| 
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Position space density. In figures [8ja)-(c) we show the time-averaged density of the 
c-field in the z = plane, i.e. 

«c(x) = <0'^(x)«Ac(x)> « ^ E '"^cCx, Qt (107) 

While the instantaneous density, e.g. see figure |7ja3)-(a4), exhibits spatial fluctua- 
tions and a random appearance of no particular symmetry, the averaged density is 
smooth and highly symmetric. The cases in figure [8];a)-(c) vary from a condensate 
fraction of < 0.5% (figure [8ja)) to 56% (figure [8]^c)), yet the spatial density profiles 
change rather gradually and do not provide clear evidence for condensation. We dis- 
cuss our procedure for quantifying the condensate later in this section. We also note 
the work of Krauth who has developed a path integral quantum Monte Carlo scheme 
for the trapped Bose gas in reference [126] and has computed density profiles for 
systems with up to lO'* atoms. 

Momentum space density. It is also desirable to be able to calculate correlation func- 
tions of the momentum space field operator, 

^c(k) = 7;rW r ^/'x^c(x)^'-*". (108) 

A particularly useful example is the momentum density, 

«c(k) = <0c(k)0c(k)>. (109) 

The use of Fourier transforms to convert from the spatial to momentum representa- 
tions of the c-field makes evaluating such observables quite efficient. 

In figure [8jd)-(f) the momentum density in the fc^ = plane is shown for the 
cases corresponding to the position space densities in figure [8ja)-(c). Noting that 
the momentum density axis is logarithmic, we observe a strong peak forming as the 
condensate fraction increases, providing an unambiguous signature of condensation. 
In figure [8jg) the momentum density for the case in figure [8jf) is shown in the = 
plane, where the anisotropy of the condensate mode (due to the anisotropy of the 
confinement potential in the xy-plane) is clearly apparent. 

Condensate. It is important to quantify the amount and nature of condensate in the 
system. Unlike the uniform system, where the condensate mode is always the zero 
momentum mode, particle interactions in the trapped system have a strong effect on 
the shape of the condensate mode and cause it to be significantly different from the 
ideal case (see equation ( [73] ) for the regime T <«c Tc). 
According to the criterion provided by Penrose and Onsager [160], the condensate 
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Figure 9. Equilibrium density for the c-field region description of a Bose gas. Condensate density (dashed black 
line) and c-field region density (solid grey line) densities along (a) jc-axis (more tightly confi ned dir ection).v and (b) 
i/-axis. Simulation for Ec = lONcTioj. with other parameters given in section [3.2.1| 

number A'^o is identified as the largest eigenvalue of the one-body density matrix, 
defined in terms of the field as 



G'^x,x') = <^c(x¥c(x')>. (110) 

The corresponding eigenvector is associated with the condensate mode of the system, 
^o(x), i.e. 

J fif^x' G'«(x,x')^o(x') = A^o^o(x). (Ill) 

For the uniform system the one-body density matrix exhibits the property of off di- 
agonal long range order (ODLRO) when a condensate is present [218], i.e., 

G'^(x,x')^y as |x-x'|^oo. (112) 

In this case the condensate orbital is ^o(x) « 1 / Vv, where V is the system volume 
and the thermodynamic Umit is assumed. 
Another widely used definition of the condensate "order parameter" is given by 

^o(x) = <(A(x)>, (113) 

based on the idea of spontaneously broken gauge symmetry. The time averaged value 
of (i/'c) in the c-field approaches is typically zero, and so this definition is of lim- 
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ited use for quantifying condensate in this review. For a comprehensive discussion 
on the various definitions of condensate we refer the reader to the review article of 
Leggett [133], who shows preference to the Penrose-Onsager definition, and states 
(on the topic of the broken symmetry definition) "... while possibly streamlining 
some calculations when judiciously used, is liable to generate pseudoproblems and 
is best avoided". 

In our formalism G'^(x, x') is equivalently and much more efliciently computed in 
the mode basis as G^f^ = {a*^a„), which is quite feasible to compute for the typical 
classical region size (< lO'* modes in C) using time-averaging. In the spectral basis 
the condensate mode is specified by a vector aj] such that ^„ Glf^a'^ = Noa'^. 

In figure|9]we show the time-averaged position density along two coordinate axes 
obtained from a c-field simulation. In addition to the total c-field density uq, we also 
show the condensate density |^o(x)P. For reference, the condensate number, A^o> for 
simulations over a wide range of energies are given in table [T] In the remainder of 
this section we develop techniques for extracting other thermodynamic quantities to 



attribute to these calculations: temperature and chemical potential in section 3.2.5 



and incoherent region atoms in section 3.2.6 



We note that the correlation functions discussed so far only apply to the c-field 



region. We return to this issue in section 3.2.6 when we consider including contri 



butions from the incoherent region. We also mention that higher order correlation 
functions, including second order (e.g. density fluctuations) correlations functions 
have been calculated using the PGPE approach, see references [19,27, 188] (also 
see [115]). A procedure for calculating two-point correlation functions is discussed 
in section U. 4. 31 

3.2.5 Thermodynamic quantities: temperature and chemical potential 

It is desirable to find a means to attribute a temperature to the thermalized state of 
a c-field simulation. Previous attempts to do this have been based on fitting the oc- 
cupation of high energy modes to perturbative calculations for the spectrum based 



on Hartree-Fock-Bogoliubov (HFB) theory [45,52] (see section 3.3.1 1. For harmon- 
ically trapped gases, calculation of the HFB modes is much more difficult, and limits 
temperature calculations to perturbative regimes. However, the temperature can be 
crudely estimated by fitting the high momentum components of the system to a non- 
interacting distribution [86]. 

An alternative approach of general applicability is found by extending Rugh's dy- 
namical definition of temperature for classical Hamiltonian systems [179] to the 
PGPE. This scheme has the advantage that it is non-perturbative, and is quite ac- 
curate. 

Rugh's approach was formulated for a classical mechanical system, and it is con- 
venient to write the c-field Hamiltonian as Hq = Hq{T), where F = {Qj,Pj} is the 
vector of the canonical position and momentum coordinates introduced in section 



2.3.6 (equations (69 1, (70 1). We also need to explicitly account for the c-field nor- 
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malization functional (62 1, Nq - Y^j |ar/| , which is another constant of motion and 
can also be written as afunction of canonical coordinates, i.e. A^c - AAc(r). The 
usual expression for the temperature of a system in the microcanonical ensemble is 
given by 

' . (1.4) 



T \dEc 
where the entropy is defined by 



S =kBln{ dT 6[Ec - HciTmNc - Nc{T)] } , (1 15) 



with the delta functions ensuring our microcanonical description is one of fixed c- 
field energy and normalization. 
There is several issues with using equation ( 114[ ) to determine the temperature in 



the PGPE approach. First, it is practically impossible to determine the entropy 5 
when a large number of modes are in the c-field region. Second, equation ( 114[ ) can- 



not be evaluated using a single micro-canonical ensemble average, since T depends 
on the derivative with respect to energy. In 1997 Rugh made a fundamental contribu- 
tion to statistical mechanics by proving that a micro-canonical average could be used 
to calculate the temperature. This approach is now extensively used in the molec- 
ular dynamics community since the micro-canonical average can be replaced by a 
time-average, as we shall do here. 
Rugh's result, proven using differential geometry methods, showed that the tem- 



perature expression ( 1 14 1 could be equivalently written as 

jL^I^O.XHd). (116) 

rigorously shown to work for Hamiltonian systems at energies where the energy 
surface is regular [179-181]. The components of the vector operator D are 

A=^'/Jr, (117) 

where e, can be chosen to be any scalar value, including zero, and the vector field 
can also be chosen freely within the constraints 

DHc-Xt = 1, DNc-'S.t = 0. (118) 



Geometrically this means that the vector field Xj- has a non-zero component trans- 
verse to the Hc(T) = Eq energy surface, and is parallel to the Nc(r) = Nq surface. 
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A vector field that satisfies these constraints is 

DHc - AnDNc 



(119) 



IDHcl^-ANiDNc-nHc) 
where we have introduced the parameter = 'DNc-'DHcl\T)Nc?'. The expectation 



value in Eq. ( 1 16 1 is over all possible states in the microcanonical ensemble and can 
be evaluated as a time-average for our ergodic c-field system. In the interests of 
brevity we do not discuss the additional technical details of how the matrix elements 



in equation ( 1 \9\ are evaluated, however point out that a procedure for doing this 



exactly and efficiently is given in reference [46]. 

It is worth giving an example to illustrate the formalism. Consider a simple system 
of M degenerate oscillators described by - Y^j^^j, with no normalization con- 
straint. Taking 2), = 5/5r, we have that (Z^), = DiHc/\DHc\^ = Ti/Hc, where we 
have used that |D//cP = Z, r] = Hq. Finally, we have that llkeT = (D ■ Xt) = 
(M - l)/Ec, which is the standard micro-canonical result. 

Similar to the discussion above, the chemical potential can be evaluated according 

to 



yu I dS 



= {D-X,iT)), (120) 



where the conditions on the vector field are 



DHc-X,^0, DNc-X,^L (121) 
The appropriate vector field is of the same form as the right hand side of equation 



( 1 19 1 but with He and A^c interchanged. 



In figure 10 we show instantaeous values of the Rugh observables for temperature 



(i.e. [kgD ■ Xt]~^) and chemical potential (i.e. ksTT) ■ X^) evaluated from a PGPE 
evolution. The time-averaged results for these parameters over a broad range of initial 
energies are given in table [T] 

3.2.6 Including the incoherent region atoms 

To relate the PGPE results back to an experimental system we need to account for the 
sparsely occupied modes of the incoherent region, which we have so far ignored. To 
do this we take the classical region and the incoherent region to be weakly-coupled 



systems in thermal and diffusive equilibrium (see figure 111, i.e. with the same tem- 
perature and chemical potential. The thermal cloud exists in the potential of the trap 
plus time-averaged c-field density ?ic(x) determined from the PGPE simulations. To 
model the incoherent region modes we use a Hartree-Fock approximation. As dis- 



cussed in section 2.3.9 the mean-field approach provides a good description of the 
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Figure 10. Extracting dynamical thermal quantities. Instantaneous value of Rugh observable for (a) temperature 
and (b) chemical potential shown as black lines evaluated over one second of evolution. Average values of T and fi 
shown as grey horizontal lines in (a) and (b) respectively. Simulation for Eq = IQNqTiu)-, with other parameters 

given in section [3.2. 1 1 




Figure 1 1 . Schematic view of the coupling between the c-field region, described by the PGPE, and the incoherent 
region. The systems are assumed to be weakly interacting and in thermal and diffusive equilibrium. 



system away from the critical region (e.g. see [82]), and of modes well-above the 
energy scale unc (also see the discussion in section 3.2.7 1. 

The average properties of the incoherent region can be calculated from the one- 
particle Wigner distribution 



Fi(x,k) = ?^ , (122) 

exp(;S[eHF(x,k) -/i]) - 1 
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where 

eHF(x, k) - — + yo(x) + 2u{nc{x) + ?ii(x)), (123) 
2m 

is the Hartree-Fock energy, and ji is the chemical potential. The one-particle Wigner 
distribution is related to the one-body density matrix for the incoherent region (see 



equation (137)), and should not be confused with the multi-mode (many -body) 



Wigner function discussed in section 2.3 



In this semiclassical description x and k are treated as continuous (commuting) 



variables. However, care needs to be taken to ensure that Eq. ( 122 1 is only applied to 
the appropriate region of phase space spanned by the incoherent region, i.e. single- 
particle modes of energy exceeding ecut- Interpreted in phase space coordinates, this 
region is 

a = |x, k : — + yo(x) > ecut I . (124) 
A quantity of particular interest for us to calculate is the incoherent region density 

r (P\i 

Jn, (2ny 



^k^Fi{x,k), (126) 
,(x) 2n'^ 



where we have made use of the isotropic nature of the kinetic energy term and have 
implemented the phase space restriction, Q.i, as a spatially dependent lower cutoff on 
the integral 

/ ^ V2m[e„,r - Vo(x)] 
^cut(x) = e{ecut - Vo(x)) , (127) 



where 6(x) is the unit step function. The incoherent region atoms interact with those 
in the c-field region, which can be accounted for by adding an effective potential SV = 
2m«i(x) to the c-field description. To lowest order this shifts the system chemical 
potential by Ayu w 2uni{0). To ensure complete self-consistency, the c-field properties 
would need to be re-simulated including the effective potential, however this is often 
unnecessary as the incoherent region density is often quite small and approximately 
uniform in the spatial region of overlap with the c-field atoms. 
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Figure 12. Total density profiles including the incoherent region, (a) Density along x-sods and (b) density along 
i/-axis. (black dashed line) condensate density, (gray solid line) c-field region density, (black dots) incoherent region 
density, (black solid line) total system density, and the pure Hartree-Fock result for the total density (grey dashed 
line). Simulation for Eq = IStioj-, with other parameters given in section [3.2.1| 

We can also calculate the momentum density of the system as 



ni(k) 



%2u2 



2m 



(128) 



Using the Hartree-Fock analysis we can now include the incoherent region atoms 
into the PGPE simulation results presented in the previous sections. In figure [T2] we 
show the typical profiles comparing the c-field and incoherent region density profiles, 
including the total density 



n{x) = ndx) + «i(x). 



(129) 



These results also allow us to ascribe the total number of atoms, = A^c +Ni, where 
Ni = J d^xni{x), to the simulated systems. Using this analysis of the incoherent 
region in table [T] we can attribute total particle number to our PGPE simulations. 
For comparison, in figure 12 we have also shown the result of a pure Hartree-Fock 
analysis (as described above, but taking ecut ^ so that all modes are treated using 
mean-field theory). The pure Hartree-Fock result is for the same temperature and 
total number as the c-field calculation, yet predicts no condensate, since the tem- 
perature is a few nK above the mean-field critical temperature. Outside such critical 
regimes the difference between mean-field and c-field calculations is generally much 
smaller. 

The results in table [T] show that for fixed c-field region (i.e. fixed cutoff tcut and 
A^c)^ the temperature and total number of particles grow rapidly as the c-field en- 
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ergy increases. In general this means to simulate a fixed total number of particles for 
various temperatures, we must appropriately manipulate the macroscopic parameters 
defining our micro-canonical system, i.e. tcut, A^c and Ec- We will see in section |5] 
that use of the stochastic PGPE formalism greatly eases the effort required to calcu- 
late systems at definite temperature by making use of a grand-canonical description. 
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Table 1. Summary of PGPE thermalizatioii results for a c-field region with Nc = 10,000 Rb-87 atoms. Other parameters: 
Ux-fij^f:} = ( 120, 30, 30( Hz and fcui - 33^w-. For reference, the Thomas-Fermi ground state energy is £tf = ^^cMtf ~ 
9.04Nctia}:,. Note: T and /j are determined by the average of two different choices of Rugh temperature (see [50]), one of which 
is shown in figure [TOj 



3.2. 7 Validity conditions 

The high mode occupancy of c-field region described by the PGPE makes the va- 
lidity requirements of this approach somewhat different from those listed for the 
truncated Wigner approach in section [2! 3. 9 In particular, the dominance of classical 



fluctuations means that the thermalization of quantum noise is not a concern. Thus, 
the conditions for the PGPE method to provide an accurate description of the c-field 
region are: 

(i) Good basis: The cutoff has to be sufficiently large that the single particle modes 
provide a good basis for describing the interacting c-field region modes. This 
condition can be expressed in terms of the peak (central) c-field density as ecut - 
eo > undO), where 6o is the ground single particle energy. This condition also 
ensures the validity of the separation into C and I regions in the critical regime 



(see discussion in section 2.3.9 1 



(ii) High mode occupation: the mean occupation of the highest energy single particle 
state is greater than unity. In general we refer to this quantity, extracted from 
simulations, as «nun which is also listed in table [T] to demonstrate the validity of 
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k (2jt/L) 

Figure 13. Fits of the simulation quasiparticle population data to dispersion relations. The dots are a plot of 
{l/Nii - i/No), the solid curve is for dispersion relation predicted by second order theory, and the dashed curve is 

the dispersion relation predicted by Bogoliubov theory. Simulation parameters: u = 200QL^£l/Nc, 
Eq = 4000A'c6i, and Nq/Nq = 0.279, where the unit of energy is = tP'/lniL?, and the unit of temperature is 
7^0 = Nc^l/^b- Reproduced from [52] © 2002 by The American Physical Society. 



those results. 



3.3 Applications to the uniform Bose gas 

3.3.1 Temperature and quasi-particle modes of the uniform system 

In a sufficiently weakly interacting Bose gas, the Hamiltonian for the system can 
be approximately diagonalized by a transformation to the Bogoliubov quasi-particle 
basis. For the uniform gas, the interaction term only mixes modes of opposite mo- 
mentum, and the transformation from single-particle modes to Bogoliubov quasi- 
particles of a well-defined momentum k depends only on the product of the inter- 
action strength u and the condensate number A^o- This is a quantity that can be de- 
termined from the PGPE calculations, and so the individual classical fields can be 
projected onto the Bogoliubov quasiparticle basis, and the time-averaged quasiparti- 
cle occupations can be accurately determined. 

When the Bogoliubov quasiparticles form a good basis, we expecte that at ther- 
mal equilibrium the c-field method will result in the mean quasi-particle occupations 



being given by the equipartition relation (100 1. If we define the condensate eigen 



value as sq = A, and require that the condensate occupation also be given by the 
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equipartition relation 



A^o = (130) 
A- fi 



then we can solve for the thermodynamic chemical potential /u = A - ksT/No^By 
substituting this result into ( |100 1 and rearranging we find 



/I 1. ^^3^^ 



ksT \Nk No)' 

where the numerator of the left hand side is the quasi-particle energy relative to the 
condensate. This suggests a prescription for determining the temperature of a simula- 
tion. The right hand side can be accurately measured by ergodic averaging in c-field 
simulations, and the left hand side can be evaluated using theoretical predictions of 
the spectrum {sk - A) with the temperature forming a single fit parameter. This pro- 
cedure was developed by Davis et al. [45,52] before the application of the method of 



Rugh for determining temperature as described above in section 3.2.5 We also note 
the non-projected classical field study of Brewczyk et al. [33]. 

Bogoliubov spectrum. In the limit of large condensate fraction Nq/Nq ~ 1, we expect 
the Bogoliubov transformation to provide an accurate description of the system, with 
the dispersion relation 



Sk- A = 



2 



+ {cfikf 



1/2 



(132) 



2m 

where c = y/Nou/mL? is the speed of sound. 

Second order spectrum. For sufficiently large interaction strengths and temperatures, 
the cubic and quartic terms of the many-body Hamiltonian that are neglected in the 
Bogoliubov transformation become important. In Ref. [143] Morgan develops a con- 
sistent extension of the Bogoliubov theory to second order that leads to a gapless 
excitation spectrum. This theory treats the cubic and quartic terms of the Hamilto- 
nian using perturbation theory in the Bogoliubov quasiparticle basis, and results in 



energy-shifts of the excitations away from the Bogoliubov predictions of Eq. ( 132 1. 

The results in figure [13] clearly show that second order theory provides a better 
description of mode occupations than Bogoliubov theory. Other results in [52] show 



' We note that the familiar result oi A = fj only holds in the thermodynamic limit, and for finite size systems it is 
often important to ensure that the chemical potential /i and condensate eigenvalue 4 are distinct. 
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Figure 14. Shift in the critical temperature of a uniform Bose gas with interaction strength determined from PGPE 
simulations with Nc = lO'" for zero scattering length. The dashed line is a linear fit to the first two data points and 
this has a slope of 1.3 ± 0.4. Reproduced from [50] © 2003 by The American Physical Society. 

that as the interaction strength increases, initially better agreement with second or- 
der theory is observed, until the validity conditions of that theory are eventually 
surpassed. In reference [50] the temperature, as determined from the spectral fitting 
procedure [52] , was shown to be in good agreement the Rugh (dynamical) tempera- 



ture (discussed in section 3.2.5 1 in the regimes where spectral fitting was valid. 



3.3.2 Shift of Tc for the uniform Bose gas 

The shift in critical temperature Tc with interaction strength for the homogeneous 
Bose gas has been the subject of numerous studies and debate for almost fifty years 
since the first calculations of Lee and Yang [131, 132]. While there is a finite shift to 
the chemical potential in mean-field (MF) theory, the shift of the critical temperature 
is zero [13]. The leading order effect is due to long- wavelength critical fluctuations 
and is inherently non-perturbative. Using effective field theory it was determined that 
the shift is 



ATc/Tco - can 



1/3 



(133) 



where n is the particle number density, a is the s-wave scattering length, and c is a 
constant of order unity [12]. Until recently results for the value of c disagreed by 
an order of magnitude and even sign, as summarised in Fig. 1 of [7]. However, two 
calculations performed using lattice Monte Carlo have settled the matter, and confirm 
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that the shift is in the positive direction with combined estimate of c w 1.31 + 0.02 
[7,119]. A number of recent improved results broadly agree, and useful discussions 
are provided by Andersen [5] and Holzmann et al. [96]. 

Here, we briefly describe the procedure used by Davis et al. [50] to calculate a 
value for c using the (uniform gas) PGPE. 

(i) For a given nonlinearity (i.e. scattering length) a randomized initial state of def- 
inite energy is evolved with the PGPE, and the temperature is determined by 



using the methods described earlier in section 3.2.5 

(ii) As the initial state energy is varied, the critical point is identified as where 
the Binder cumulant, Cb = {Nq)/{No)^, with A^o the population of the zero- 
momentum condensate mode. This Binder cumulant characterizes condensate 
number fluctuations, and takes the universal value of C™' = 1.2430 at the transi- 
tion. 

(iii) The shift in the critical temperature is calculated as a function of interaction 
strength, parameterized by the s-wave scattering length. 



By fitting a straight line to the first two points as illustrated in figure 14 we get an 
estimate for the coefficient 

c- 1.3 + 0.4, (134) 

where the error specified is due to the uncertainty in the value of Tc for the data point. 
This agrees with the value determined in Refs. [7, 1 19]. 



3.4 Applications to the trapped Bose gas 

3.4.1 Shift in Tcfor a trapped Bose gas: Comparison with experiment 

The behaviour of Tc for the harmonically confined Bose gas is drastically different 
from the uniform gas. There is a shift in Tc due to finite size effects arising from 
the fact that the system is not in the thermodynamic limit [89], and a first-order 
interaction shift due to mean-field effects [84]. 

For a typical BEC experiment, the critical temperature deviates from the ideal gas 
result only by a few percent. Thermometry of Bose gases at this level of accuracy is 
challenging: However, in 2004 Gerbier et al. reported precise measurements of the 
critical temperature for a range of atom numbers [82]. 

Davis et al. [47] used those measurements to make the first quantitative compari- 
son of the PGPE formalism with experiment and other theories, which are summa- 
rized in figure 15 The various other theories appearing in figure 15 are: 

Al: This is the meanfield analytic estimate as calculated by Giorgini et al. [84], 
and was compared with the experimental data in [81]. 

A2: This is the analytic estimate as calculated by Arnold et al. [8], which includes 
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Figure 15. Compaiison of theoretical calculations with experiment. The main figure plots Tc vs Nq, whereas the 
inset plots the shift of Tc against the relevant small parameter a/zio- Experimental results: data (open circles), one fit 

(gray area). Theoretical results for T^: ideal gas (dot-dashed line), Al (dotted line), A2 (dashed line), MF-GPE 
(crosses), MF-HFBP (dots), PGPE (pluses). Solid lines through the data points are polynomial fits. Al is not shown 
in the main figure for clarity. The total number of atoms at the critical point is N = 4.0 X lO'' and 
Aq = hi ^llmnkisT. fits to the data. Reproduced from [47] © 2006 by The American Physical Society. 



next order fluctuation results, however is only strictly valid in broad traps. 

MF-GPE: The GPE is solved numerically using a variational Gaussian ansatz, and 
the thermal cloud calculated using a semi-classical approximation [84] . At each tem- 
perature the condensate and non-condensate are determined self-consistently with a 
fixed number of particles, and the critical temperature is taken to be where the con- 
densate fraction decreases to zero. This approach difl'ers from theory Al because it 
avoids using perturbation theory around the ideal (saturated) gas density profile to 
estimate interaction effects. 

MF-HFBP: Here the condensate fraction is fixed, and the temperature determined 
that gives an appropriate self-consistent condensate mode and thermal density (The 
full Bogoliubov modes are used and the semiclassical approximation is avoided). We 
have verified the results are unchanged for equipartition or Bose -Einstein statistics. 

The PGPE calculations appear to provide the best theoretical description of exper- 
iment, however, error bars in the experimental results are not yet small enough to 
definitively discriminate between results. 

We also note that the Hartree-Fock-Bogoliubov Popov calculations (MF-HFBP) 
use the same procedure as in PGPE calculation to determine the critical point, the 
above cutoff density and the total atom number, so that the difference between this 
best meanfield calculation and the PGPE results is due to beyond meanfield fluctu- 
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ation effects. This suggests if experimental accuracy in thermometry could improve 
by an order of magnitude, then effects of fluctuations on the critical temperature in 
this system could be directly investigated. 

3.4.2 Quasi-2D Base gas 

The phenomena of superconductivity and superfluidity are striking manifestations of 
the role played by quantum statistics at low temperatures. Altering the temperature 
or effective dimensionality may radically change the physical properties of quantum 
degenerate systems. A well known consequence is that in contrast to 3D, there is no 
BEC for a homogeneous 2D ideal-gas in the thermodynamic limit at any finite tem- 
perature [94, 140]. Nevertheless, the Berezinskii-Kosterlitz-Thouless (BKT) vortex 
binding-unbinding phase transition allows the emergence of superfluidity in 2D sys- 
tems [14, 124]. Although weak particle interactions alone are not sufficient to change 
the situation, an external confinement modifies the density of states in such a manner 
that the critical point of BEC is elevated to a finite temperature [9]. Therefore it is 
not certain a priori whether the transformation from normal to superfluid in such 
systems is a BEC or BKT-type transition. 
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Figure 16. Interference pattern (a) produced by two independent c-fields (b) and (c) at temperature T = 0.86 Tq. 
The relevant particle numbers are N^i = 3.0 X lO' and N = 4.0 X lO''. The "zipper" structure in (a) is the telltale 
signature of the phase singularity associated with the central vortex in (b). The locations of vortices and antivortices 
are marked by + and - signs, respectively. Reproduced from [190] © 2006 by The American Physical Society. 



This properties of the finite temperature trapped 2D system have proven difficult to 
analyze. Strong fluctuations mean that meanfield approaches are inapplicable, how- 
ever since these fluctuations are classical in nature the PGPE approach is appropri- 
ate. Indeed, early studies of the uniform 2D Bose gas were performed by a c-field 
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method, but sampled using Monte Carlo techniques [172, 173]. More recently Gior- 
getti et al. [83] have developed an accurate semi-classical approach for simulating 
this system. 

Simula et al. [190] used PGPE simulations of quasi-2D Bose-fields to characterize 
the low temperature phases for such systems over a wide parameter range. These 
simulations show the emergence of thermally activated vortices (see figure [T6jb)- 



(c)), their influence on interference patterns (see figure 16 a)), and provide strong 
evidence supporting the view that the BKT-type phase was observed in the recent 
experiment by Stock et al. [200]. More recent work with PGPE simulations have 
characterized correlation and collective mode properties of the system to quantify 
the BKT transition point [21, 188], and find good qualitative agreement with recent 
quantum Monte Carlo simulations [97]. 

3.4.3 Two point correlation functions 

Recent experimental developments in ultra-cold gases [70,88, 109, 157, 178, 183] have 
allowed atomic correlation measurements that are analogous to the photon correla- 
tions observed in the landmark experiments of Hanbury-Brown and Twiss [91]. Such 
correlations are of particular interest in systems where many-body interactions are 
important [4, 208] , and in the region of the phase transition, where critical exponents 
can be measured [59]. 

The PGPE description is valid in this regime and can be used to calculate these 
correlations, and assess beyond meanfield effects (c.f [57, 146]). We now summarize 
an approach for calculating these correlations within the PGPE formalism that has 
been developed by Bezett et al. [19]. 

The quantities of interest are the normally ordered first-order correlation func- 
tion, G^''(x, x') = (i/''(x)0'(x')) (also known as the one-body density matrix), and 
second order correlation function, G^^'(x, x') = (i^'(x)i^'(x')0'(x')i/'(x)). From these 
functions first- and second-order observables can be directly obtained, such as the 
density-density correlation function. 

Breaking the full quantum field into c-field and incoherent parts, the correlation 
functions can be written as 

G<'\x,x') - G^'^(x,x') + gJ'^(x,x'), (135) 
G<2\x,x') - G^^'(x,x') + Gf (x,x') + 2G['\x,x')G^''(x,x') 

-Hni(x)?ic(x') + ni(x')«c(x), (136) 

where G;.'\x,x') = (.a;(x')«A/x)) and Gf (x,x') = (0';(x')«A;(x)«A./(x)(A/x')) with 
i = {I, C} for the incoherent and classical regions respectively, and we have neglected 
any correlations between the c-field and incoherent regions. 
While the c-field correlations can be evaluated using the approach detailed in sec- 



tion 3.2.4 for the incoherent region we can make use of the one-particle Wigner 
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Figure 17. (color online) Two point position space correlation functions of a hamionically trapped Bose gas of 
W = 3 X lO' *''Rb atoms atT = \59nK. Other parameters: Nq = 3540, = i^Hwy, with 
{6J.V, c^y, tJ;) = 2;r { 1, 1, VS) X 40 i"' . Reproduced from [19] © 2008 by The American Physical Society. 



function given in equation ( 122| ). Appropriately transforming the Wigner function 



we obtain the first order correlation function, i.e. 

gJVx')-J d\e-'^-(^-^'^ Fii^^^,kj. (137) 

As the Fi Wigner description of the incoherent region is Gaussian, we can easily 
obtain the second order correlation function g'^\x,x') = ni{x)ni{x') + \g[^\x,x')\^. 

Figure [TT] illustrates correlation functions for a trapped Bose gas at J ~ Tc. The 
results shown are for the case of two points along the x-axis of the system, and in 
figure p^b) and (c) the normalized correlation functions, defined as g''^\x,x') = 
G'-^\x,x')/ y/n{x)n(x') and g'^^\x,x') = G^^\x,x')ln{x)n{x'), are shown. The broad 



feature apparent in figure 17 a) and (b) is the off diagonal long range order, arising 
from the emerging condensate in this system. The diagonal ridge is due to short range 
thermal correlations. We also note that Holzmann et al. [95] have used a quantum 
Monte Carlo method to obtain the pair distribution function for a trapped Bose gas. 



3.5 Applications of non-projected classical fields at finite temperature 

As well as the work of the current authors on quantitative projected c-field tech- 
niques, there have been a number of other studies of finite temperature properties of 
degenerate Bose gases. For completeness, in the final part of this section we briefly 
describe the results that have been obtained. 

3.5.1 Homogenous gas 

In one of the first papers on classical fields, Goral et al. [85] demonstrated the ther- 
malization of a homogeneous multimode Bose gas in a similar manner to Davis et 
al. [51,52]. They expanded the equation of motion for the c-field in terms of mode 
coefficients, and calculated the nonlinear terms by performing the appropriate sum- 
mations, and implicitly correctly applied a projection operation. Brewczyk et al. [33] 
performed a Bogoliubov analysis of homogeneous Bose gas using a GPE classical 
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field description. Zawitkowski et al. [220] attempt to describe not just the classical 
modes but the I region modes using only the GPE by fixing a grid cutofi" such that 
the condensate fraction and temperature agree with that for the ideal Bose gas. Doing 
this eliminates any possibility of describing e.g. describing the effects of interactions 
on the transition temperature, and includes a large number of modes in the problem 
that should not be described classically. It should be clear from this review that the 
current authors strongly disagree with this approach. 

Leadbeater et al. [129] studied the effect of condensate depletion on the critical 
velocity when an object is dragged through a superfluid. On a related note, Zaw- 
itkowski et al. [221] performed an interesting study of placing a homogenous moving 
condensate in a static thermal cloud, and investigated the decay of the superflow as a 
function of velocity and temperature. Unfortunately it seems that the lack of projec- 
tion caused some numerical issues in this work, such as the violation of momentum 
conservation. 

Witkowska et al. [216] related the dynamics of a nonlinear string to the weakly 
interacting Bose gas. Nunnenkamp et al. [150] made a comparison of three versions 
of a classical field theory for a one-dimenional Bose gas on a ring. They found that 
an exact solution in the high temperature limit of a transfer integral method agreed 
well with both a molecular dynamics approach, and classical field simulations of 
the GPE. Recently Sinatra et al. [192] found nondiffusive phase spreading of a 3D 
homogenous Bose-Einstein condensate at finite temperature. 

Connaughton and co-workers have studied condensate formation in the homoge- 
neous gas using a GPE model [41, 1 13]. In an interesting application related to classi- 
cal fields, Picozzi and co-workers have investigated the dynamics of equihbration in 
incoherent nonlinear optics both theoretically and experimentally. See, for example, 
references [127, 163-165]. 



3.5.2 Trapped gas 

Goral et al. [86] were the first to apply the condensation criterion of Penrose and 
Onsager [160] to a classical field. They solved a non-projected GPE for the trapped 
Bose gas at finite temperature, and developed some estimates of thermodynamics 
properties of the system. Schmidt et al. [184] applied the same simulation technique 
to investigate the decay of an off-centre vortex in a harmonic trap at finite tempera- 
ture. Recently Gawryluk et al. [80] seeded a trapped F = 1 spinor condensate with 
thermal fluctuations and studied the resulting spin dynamics. 

The effect of thermal fluctuations in Bose gases is more significant in low dimen- 
sions. Kadio et al. [115] studied the coherence properties in a quasi- ID trapped Bose 
gas, and analysed the effects of phase fluctuations in a 3D elongated trap. Mebrahtu et 
al. [139] have analyzed coherence effects in the spatial splitting of a quasi- ID Bose- 
Einstein condensate and its subsequent merging at finite temperature. 
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3.5.3 Superfluid turbulence 

Finally, we mention work in the area of superfluid turbulence, which attempts to 
describe the formation of tangles of vortices in the homogeneous superfluid transition 
and the subsequent relaxation to global phase coherence. Many of the simulations of 
these systems make use of the GPE to describe finite temperature non-equiUbrium 
dynamics, and hence are directly related to classical field techniques. 

Of particular interest to this review is the work of Berlofl" and Svistunov [16], 
who studied condensation from a strongly nonequilibrium state in a 3D homoge- 
neous system using the GPE. Their main interest was in the decay of superfluid 
turbulence, and the establishment of phase coherence, validating the scenario of su- 
perfluid growth as earlier described by Svistunov and Kagan [116-118,205]. Berloff 
subsequently studied the interactions of vortices and solitary waves and their role in 
the decay of superfluid turbulence [15], as well as turbulence in a two-component 
system [17]. Recently Berloff and Youd studied the decay of vortex rings in a ho- 
mogeneous superfluid at finite temperature [18]. In reference [122] Kobayashi and 
Tsubota simulate a GPE with a dissipation at short wavelengths and obtain an energy 
spectrum consistent with the Kolmogorov law. 



4 Applications of the truncated Wigner PGPE to quantum matter wave dynamics 

As more experimental investigations have begun probing beyond mean field quantum 
dynamics in BECs, theoretical applications have begun to explore the role of thermal 
and quantum fluctuations using the truncated Wigner method. Here we give a brief 
survey of the background and recent developments of this method. 

Background. The truncated Wigner approximation was introduced by Robert Graham 
in 1973 [87] and has met with wide application in the field of laser physics. A precur- 
sor of work on trapped Bose gases was carried out by Carter et al. [36] who applied 
phase space methods to the simulation of the quantum optical nonlinear Schrodinger 
equation. The theoretical formulation and applications to Bose gas dynamics began 
with the work of Steel et al. [197] who developed phase space techniques for atomic 
Bose fields and applied them to simulating the time evolution of a one dimensional 
homogeneous Bose gas. The truncted Wigner method was compared with the func- 
tional positive-P phase space method [61] in calculations of the first order coherence 
function g^^\t) = (aJ(Oao(0))/(aJ(0)ao(0)> for the condensate operator ao- Different 
initial states of the condensate were sampled including the coherent state and the Bo- 
goliubov state. A general conclusion of this work, which provides a reliable guide, is 
that the positive-P method, while exact, is unstable except for very short simulation 
times, whereas the truncated Wigner method, while approximate, is stable. Many 
subsequent works have considered aspects of the validity of the truncated Wigner 
method and its applications to dynamical Bose gases using both full phase space 
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Figure 18. Schematic of condensate collision scenario of [147]. (a) Position space densities of initially 
overiapping, counter-propagating condensate wavepackets. (b) Momentum space representation of possible energy 
and momentum conserving collisions between atoms in the two condensates (kj, k2) onto the allowed spherical 
scattering halo (ks, k4), indicated by the gray annular region in the collision plane. 

approaches, and the classical field method based on analysis of single trajectories. 

A distinction between the formulation of [197] and the c-field description pre- 
sented in section[2| is the emphasis placed on projection into a low energy subspace, 
both formally and numerically. In the TWPGPE formulation the projection operator 
imposes a formal UV-cutoff which allows a measure of control over the sometimes 
spurious effects of vacuum noise arising in the truncated Wigner method. 



4.1 Condensate collisions in free space 

The Bragg scattering of a condensate into a superposition of states of momentum 
and 2fik creates a well characterized non-equilibirum initial condition that is easily 
produced in experiments [25, 125, 198, 199]. This scenario is shown schematically in 



the centre-of-mass (COM) frame (in position space) in figure 1 8 1, where the original 
and scattered wavepackets move away from each other. In the subsequent dynamics, 
but while the two wave packets still overlap in position space, pairs of atoms are 



scattered onto a spherical shell in momentum space (see figure 18 d). This scattering, 
often referred to as an S-wave halo, is clearly seen in experiments [39, 161], but is 
absent in a GPE description. 

The truncated Wigner method was first used to model this process by Norrie et 
al. [147, 148]. Beginning with a condensate with mode function ^o(x), Bragg scatter- 
ing was assumed to scatter half the condensate, resulting in the superposition of two 
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Figure 19. (Color online) (a)-(f) Velocity mode populations on the planes v- = (left) and Vj^ = (right) for the 
condensate collision described in the text at ; = (top), t = 0.5 ms, and t = 2.0 ms (bottom). The spherical 
momentum cutoff is clearly visible in the upper plots due to the presence of quantum flucuations. (g)-(h) Mode 
populations at f = 2.0 ms for an identical collision excluding vacuum noise. Reproduced from [147] © 2005 by The 

American Physical Society. 



wavepackets with momenta +/zk (in the COM frame) [24], i.e., 



(138) 



The full initial condition (see figure 19i-b) was sampled by adding vacuum noise 



to modes orthogonal to i^o (see equation (90 1). In the truncated Wigner simulation 
modes on a spherical shell of radius i; ~ 10 mm/s in velocity space are seen to grow. 



while the initial wavepackets are situated at the poles of this sphere (see figure 19:- 
d). The importance of vacuum fluctuations is clear: they seed the growth of the halo 
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Figure 20. Reduction in condensate population during a condensate collision. Wigner simulations (solid), a linear 
fit to the rate given by Fermis golden rule (dotted) and the solution to the differential equation equation \139\ that 
includes condensate depletion (dashed). Reproduced from [68] © 2008 by The American Physical Society. 



modes, thus mimicking spontaneous processes. It was found that after the halo first 



develops (see figure 19 >f), the stimulated evolution of these scattered modes leads 



to turbulent dynamics. In contrast, no such spherical shell is seen to develop in the 
GPE simulations (see figure 19 j-h). 

We also note the work of Deuar et al. [55] which used a positive-P method 
(see [42, 54, 60]) to simulate BEC collisions. The positive-P method is a phase space 
approach, like the Wigner method, but does not require approximation (i.e. the trun- 
cation of 3rd order derivates) to arrive at a set of stochastic equations for any Hamil- 
tonian which is at most quartic in operators. The application to condensate collisions 
represents a significant success of the method for modeling real time dynamics of 
atomic Bose fields. While providing an exact mapping and being used widely for 
treating quantum optical systems (where dissipation is usually significant and inter- 
actions weaker), for pure Hamiltonian evolution the method suffers from stability 
problems limiting its use to short simulations and low density systems. Comparing 
with truncated Wigner simulations, Deuar et al. found discrepancies between the two 
approaches for the early time dynamics of initially unoccupied modes into which 
atoms were scattered. While directly probing such a discrepancy in an experiment 
would be difficult, this result emphasizes the need for care in interpreting TWPGPE 
results for lowly occupied modes. 



Condensate depletion. A comparison between the truncated Wigner method and 
Fermi's second golden rule was made by Ferris et al. [68] for the case of colliding 
condensates in the uniform system. That study compared the early time depletion of 
the colliding condensates due to spontaneous scattering with the Fermi golden rule 
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prediction 



dt ~ InffiV ^^^^ 



where A^o is the number of remaining (unsc altered) condensate atoms, and V is the 
system volume. The condensate population from the truncated Wigner simulation is 



shown in figure 20 where the good agreement with the Fermi golden rule estimates 
is evident at short times. The growing discrepancy at long times arises from the 
depletion of the condensate and the stimulated dynamics of the scattered modes. 



4.2 Truncated Wigner treatment of three-body loss 

The three-body loss process is an inherently non-diffusive process in phase space 
(the generalized Fokker-Planck equation contains derivatives beyond second order 
in the phase space variables), and thus does not admit an exact formulation in terms 
of stochastic differential equations for any pseudo-probability distribution. 

An important consideration must be borne in mind at this point: practical applica- 
tion of stochastic methods is not feasible for problems containing higher than second 
order derivatives for phase space variables. Thus, a more general statement of the 
TWA is that it should include all terms up to second order in phase-space variable 
derivatives, subject to the usual validity conditions of weak interactions and high 
mode occupation. This is the approach taken in the treatment of three-body loss. 



The basic technical extension beyond the standard TWA presented in section 2.3.6 
is an additional stochastic term which models the diffusive effects of inelastic loss. 
Thus, while elastic two body collisions do not generate a stochastic equation of mo- 
tion within TWA (second order terms ai^e identically zero), three body inelastic col- 
lisions introduce a stochastic element to the evolution. 

The three-body loss master equation for the system density operator p, which has 
been rigorously derived by Jack [103, 104], takes the form 



dp 
dt 



^ j d\ {2^ix)'p{fr\x)' - ^+(x)3,A(x)'p - p^\xf4r{xf} , (140) 
which generates the time evolution for the total atom number 

^ = -Ki j d'x g3(x)n(x)\ (141) 



where 



(.A'(x)^0'(x)') 
^3(x) = . . ■ (142) 
<(A'(x)iA(x))3 
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Within the truncated Wigner approximation equation ( 140 1 leads to a stochastic dif- 
ferential equation for the c-field i/^cCx) 



#c(x) = Tc |-y |^c(x)lVc(x¥f + y ^|(Ac(x)rfi?W3(x, 0| , (143) 

[149] (in addition to the terms already in the TWPGPE ( [60] )) where the noise term 
is given b}Q 

dW^{x,t) = Y,d^n'Pn{^), (144) 

neC 

with d^n{t) a complex Gaussian noise satisfying 



dUmn'it) = 0, (145) 



d^;(t)d^n'{t) = 5„„,dt. (146) 

Application to condensate collapse. The three-body loss formalism was originally de- 
rived and applied to quantify the atom loses in condensate collisions [147, 148], 
where it was shown to be a small effect. A regime where three-body corrections 
are more important is in the description of the Bose-nova experiment performed by 
Donley et al. [58]. In that experiment a Feshbach resonance was used to suddenly 
change the scattering length from a value of a « (ideal stable BEC) to a nega- 
tive value (i.e., attractive interactions), causing the system to collapse. During this 



process the condensate density increases significantly, and from equation ( 141 1, it is 
clear the three-body loss will become more important. This problem was studied with 
the TWPGPE approach by Wiister et al. [217], who assessed the effects of quantum 
and thermal fluctuations on the collapse process. Where comparison was possible, 
the TWPGPE simulations of the collapse process agreed quantitatively with the re- 
sults of Hartree-Fock-Bogoliubov theory, and both theories predicted slower collapse 
than observed in the experiment. 



4.3 Quantum reflection of a Bose-Einstein condensate 

Scott et al. [185, 186] used the GPE and the truncated Wigner approximation to 
model the collision of a BEC with an abrupt potential barrier, as studied experimen- 
tally by Pasquini et al. [158, 159]. The system consists of a BEC held in a magnetic 



Note that while this SDE gives an exact con'espondence with the TWA Fokker-Planck equation derived from equa- 
tion )140) , the noise must be computed on a dense quadrature grid to ensure that the SDE preserves the correspon- 
dence (see Appendixlclfor a discussion of quadrature methods and SDEs). 
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Figure 21. Experimental absorption image of BEC for an impact velocity of = 3.0 mm at f = 120 ms, 
having reflected from the Casimir-Polder potential of a pillared silicon surface. The field of view is 500 fim, the 
vertical dashed line indicates the position of the barrier Lower inset: corresponding simulated absorption image in 

the y - X plane including quantum fluctuations for reflection from a barrier of height V = 1.67 X 10"" J. Upper 
inset: equivalent constant density surface excluding quantum fluctuations, axes are shown in the figure. The BEC in 
this simulation has a peak density of 5.2 X lO'-cm"' with its long axis perpendicular to barrier. Reproduced 
from [185] © 2006 by The American Physical Society. 



trap which is then accelerated at normal incidence toward a steep potential drop. Two 
regimes of behaviour for these reflections were characterized: (i) For low approach 
velocities the BEC was observed to suffer disruption due to the interference of inci- 
dent and reflected components. Most aspects of these slow collisions were adequately 
explained by the GPE, however for dense initial condensates the inclusion of vacuum 
fluctuations was observed to have an appreciable efi'ect on the dynamics through the 
formation of a scattering halo (see figure [2T]). (ii) At higher velocities there is neg- 
ligible disruption due to interference, so that the GPE results are relatively smooth. 
Studying this regime with the truncated Wigner approach, the inclusion of vacuum 
fluctuations cause a large scattering halo to develop. 

In both regimes the experiments and the truncated Wigner results were found to be 
in quantitative agreement. 



4.4 Applications to optical lattices 

There has been several studies of atom dynamics in ID optical lattices using the 
truncated Wigner approach. 
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Dynamical instability of a BEC at the band edge of an optical lattice. Ferris et al. [68] 
presented experimental results and truncated Wigner simulations of the dynamically 
unstable evolution of a BEC prepared in a band-edge state in a ID optical lattice. 

The theoretical description was based on a full 3D simulation of the experimental 
system, and included degrees of freedom transverse to the lattice, and excited band 
states along the lattice direction. The large number of modes needed to accurately 
model the actual experimental system (i.e., in the combined lattice and weak har- 



monic potential), would violate the validity condition A^c ^ (see section 2.3.9 1. 
To avoid this, the theoretical model was simplified to a translationally invariant case, 
greatly reducing the number of basis modes required. The truncated Wigner simula- 
tions showed that vacuum fluctuations have an important role in seeding the growth 
of unstable modes, leading to rapid depletion and heating of the condensate. Further- 
more, the drastic modifications of energy and momentum conservation in the lattice 
were observed to have a substantial effect on the initial dynamics in the system, par- 
ticularly the modes into which atoms were spontaneously scattered. 

Quantum fluctuation effects on dipolar oscillations. In [169] Polkovkinov et al. con- 
sidered the dipolar motion of a condensate displaced relative to the centre of the 
harmonic trap in a quasi- ID lattice. This study, conducted within the tight-binding 
Bose Hubbard description [26, 108], examined the nature of the damping, and how 
it is influenced by the quantum fluctuations included within the truncated Wigner 
treatment. An adiabatic mapping procedure was used to sample the initial Wigner 
distribution which allowed them to sample the ideal Bose gas in the harmonic and 



lattice potential [29] as discussed in section 2.3.8 Their study presented evidence that 



there is a smooth crossover between the classical localization transition (overdamped 
oscillations that occur beyond a critical displacement) and the superfluid-to- insulator 
quantum phase transition in the limit of zero trap displacement. Using a similar tight 
binding approach, the evolution of phase coherence in a deep quasi- ID lattice with a 
large number of atoms per lattice site was examined and compared with experiments 
by Tuchman et al. [210]. 

Number squeezing in ID lattices. Ruostekoski and coworkers have also considered 
quasi- ID optical lattice systems using the truncated Wigner approach, but used a 
beyond tight binding description [102, 182], which included excited band states. 

In reference [102] they consider the effect that lattice loading has on a quasi-lD gas 
initially prepared in a harmonic trap. The initial state was sampled using the Bogoli- 



ubov procedure [177,214] (see section 2.3.7), and then evolved through a simulated 



lattice loading procedure. Coherence and number fluctuations were evaluated, and 
observed to be in qualitative agreement with the experiments of Orzel et al. [156]. 

In later work Ruostekoski et al. considered the quantum dynamics in shallow lat- 
tices [182], and modelled experiments by Fertig et al. [69] of dipolar motion of a 
BEC in an optical lattice. In their simulations the initial state was sampled using the 
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Bogoliubov procedure for the combined harmonic and lattice potential. Using the 
truncated Wigner approach they modeled the sudden trap displacement and subse- 
quent dynamics and found qualitatively good agreement with the damping behaviour 
observed in experiments. These results, which are for the low atom number regime 
(i.e. A^c ~ 10^)> where the strict validity conditions (see section 2.3.9 1 for the Wigner 
approach are not satisfied, provides an indication that the Wigner method has an 
extended range of applicability. 



Dephasing in ID interferometers. Bistritzer et al. [22] examined the effect of quantum- 
phase fluctuations on a condensate split into two parts. Their system consisted of a 
pair of quasi- ID Rose gases which realize a basic model for an atom interferometer. 
These ID gases were modeled with a Bose Hubbard Hamiltonian using the Wigner 
method (note there was no explicit optical lattice potential, but the gases were treated 
in a lattice approximation). In detail they used the adiabatic mapping procedure to 
sample the initial Wigner distribution as discussed in section 2.3.8 In this application 
an initial non-equilibirum state was prepared by imposing a relative phase between 
the two quasi- ID systems. The dephasing between the (uncoupled) systems was then 
calculated as a function of time and shown to decay exponentially with little sensi- 
tivity to temperatures below a characteristic temperature T* . 



Quantum phase transition in a ID optical lattice. Dziarmaga et al. [65] investigated the 
quantum phase transition from Mott insulator to superfluid transition in a periodic ID 
optical lattice. The Kibble -Zurek mechanism (KZM) predicts that when the lattice is 
ramped down in a time of order tq (ramping up the tunneling rate) the winding num- 
ber will grow through a random walk of EEC phases. In the high occupation regime 
the TWA leads to a discrete nonlinear Schrodinger equation of motion, and the initial 
conditions were sampled to approximate a Mott insulator state \n, n,n, . . ., n), where 
n is the (integer) number of particles per lattice site. The KZM scaling of the winding 
number with tq was confirmed using truncated Wigner simulations of the transition 
dynamics. 



4.5 Dynamical instabilities and quasiparticle dynamics 

Quantum de Laval nozzle. The quantum dynamics of a dynamically unstable super- 
sonic current have been investigated by Jain et al. [106]. They considered the sta- 
tionary flow of a condensate in a ID toroidal trap that was modified to form a double 
de Laval nozzle geometry by the inclusion of a spatially varying potential (around a 
torus of length L) of the form 



V{x) = - 




(147) 
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Figure 22. Stationary flow in a doubly constricted toroid potential: (a) Schematic of a hydrodynamic de Laval 
nozzle. A flow that attains the speed of sound (v = c) at the narrowest point of the nozzle becomes supersonic 
beyond the nozzle, (b) Stationary GP solution in a quasi-one dimensional toroidal geometry perturbed by the 
potential {l47j, with winding number wq = 10 (other parameters are given in [106]). The full solution is compared 

with the results of hydrodynamic and perturbation theory approaches showing the importance of beyond 
hydrodynamic corrections in the supersonic region, (c) Flow velocity and sound velocity corresponding to the GP 
solution of (b). showing the position of the black hole (BH) and white hole (WH) horizons for sound waves. The 
energy unit of the toroid is Sto/, = ti^/mL^. Reproduced from [106] © 2007 by The American Physical Society. 



which has periodicity 2 over the region -L/2 < x < L/2. There are persistent current 
solutions for this system which have distinct spatial regions of subsonic and super- 
sonic flow, with two acoustic horizons for sound waves (phonons): one where the 
flow goes supersonic (black hole, see figure 22 1) and the other where it returns sub- 
sonic (white hole). A typical stationary solution of the GP equation which has this 
property is shown in figure [22j5 and the flow scenario is shown in figure [22]: . Super- 
sonic flows are known to be energetically unstable and will decay in the presence of 
dissipation unless the decay is topologically prohibited such as in a toroidal trapping 
configuration. It was found that the system can also be dynamically unstable in cer- 
tain scenarios, and the quantum dynamics of the instability were investigated with 
the truncated Wigner approach and compared with the predictions of Bogoliugov 
theory. 
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Figure 23. Mode populations from averaging 40 trajectories of the truncated Wigner evolutions for the Quantum 
de Laval nozzle. Modes 1 and 6 are dynamically unstable, corresponding to the negative and positive energy modes 
respectively. Reproduced from [106] © 2007 by The American Physical Society. 

For a system with a dynamical instability the usual Bogoliubov expansion (72) is 
inadequate [106, 134], as modes arise with complex eigenvalues. For these unstable 
modes the Bogoliubov description acquires an irreducible off-diagonal component, 
which takes the form (equation (47) of [106]) 

j ' ' ^ 

-2 / nyj[{bjibj.-b]^b]j + J dxiUjVj - Uj*Vj*)] (148) 



+ 



where [UJ, Vj] and [U ^ , V ■ } are the positive and negative energy Bogoliubov modes 
respectively, dr = (ioij - \y-^fi are the respective (complex) eigenvalues for these 

modes with annihilation operators [134]. The second line is analogous to the 
interaction Hamiltonian for non-degenerate down-conversion of light by a nonlinear 
crystal. In regimes where only a pair of modes are coupled then Hi describes the 
formation of a two-mode squeezed state [21 1]. In this case, tracing over the negative 
energy mode, the density matrix for the positive energy state is of the form of that 
for a thermal state with mean occupation («)+ - sinh^ yt. This is analogous to the 
Hawking effect in that pairs of quasi-particles are produced at no energy cost: one 
enters the negative energy state (in the supersonic regime) and the other is promoted 
to positive energy and emerges at the horizon. 

While the Bogoliubov analysis is useful in predicting the regions of instability, it 
cannot describe the process dynamics, as the unstable modes grow exponentially (at 
least initially) and rapidly invalidate the linearized Bogoliubov analysis. Simulations 
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Figure 24. (Color online)(a) Classical field density of equilibrium state of stin'ed condensate. Vortices are 
indicated by + symbols, and are indicated only where the surrounding density of the fluid exceeds some threshold 
value, (b) Power spectral density traces at particular radii. The black (blue) and dark grey (red) lines con'espond to 
radii r = 3.1894ro and r = 11.9575ro respectively. Data coiTesponds to the period t = 9900 - 9910 trap cycles. The 
plots in (a) and (b) are from the same simulation, in which the condensate was initially in a trap of frequency oj,-, at 
temperature T = 0, chemical potential /i, = XAfiLj,-. The elliptical perturbation rotates continuously at n = 0.75(ii,.. 

VS/njaif. Reproduced from [207] © 2008 by The American Physical 
Society. 



The spatial scale is the oscillator length ro : 



with the truncated Wigner approach (see figure 23 1 avoid such limitations as they 
include the nonlinear interactions between excitations and their back-reaction. These 
features of the Wigner approach have seen several recent applications to the study 
of cosmological analogue models in BEC systems, such as particle production in an 
expanding universe [107], and studies of Hawking radiation [37]. 

In relation to the treatment of instabilities, also note the work of Polkovnikov [170] 
on the evolution of the macroscopically entangled states in optical lattices, where the 
truncated Wigner approach was used to deal with an unstable system where the usual 
Bogoliubov treatment breaks down. In that work the author presents arguments that 
this formalism should be able to adequately describe collapse and revival dynamics 
of the condensate. 

Another application to quasiparticle dynamics was performed by Modugno et 
al. [142] who investigated the possibility of driving a parametric resonance in a 
toroidally trapped BEC. It was shown that specific quasiparticle modes could be 
resonantly excited, and then individually observed via expansion imaging. Starting 
from a zero temperature BEC, the quasiparticle excitation was initiated from vac- 
uum fluctuations present in the initial state of the Wigner representation and driven 
by modulating the trapping potential. 
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4.6 Vortex formation in a stirred BEC 

A number of experiments have shown that rotationally stirring a low temperature 
BEC can lead to formation of a vortex lattice e.g. [1,2,93, 137]. In the most typical 
scenario, the stirring excites a dynamical instability of the condensate [195], which is 
transformed into a highly turbulent state, and then after a long period evolves into a 
rotating state containing a regular vortex lattice. This system provides a challenging 
test for dynamical theories of cold bosonic gases, because while the process (ideally) 
involves only an initial pure ground state subject to a conserving (Hamiltonian) pro- 
cess, dissipation is required for the system to evolve into a state in equilibrium with 
the stirrer, e.g. see [114,209]. A number of approaches have been given based either 
on the pure GPE, e.g. [136, 145] or the GPE supplemented with phenomenological 
damping terms [209] , but the description of the turbulent state is clearly beyond the 
validity of the GPE. An a priori description of the thermalisation which is central 
to the process, and provides the necessary dissipation mechanism, was first given by 
Lobo et al. [135] using a classical field method. They simulated a three dimensional 
condensate stirred by an elliptically perturbed rotating harmonic trap, and considered 
the case of an initial T - condensate, (which they modelled as the ground state of 
the GP equation), and also the case of an initial finite temperature condensate. Their 
simulations for an initially T - Q condensate showed evolution similar to that seen 
in earlier approaches (e.g. [1 14, 136, 145]), with vortices eventually entering the high 
density region of the field a few hundred trap cycles after the creation of the turbulent 
state. The vortices settle into an ordered lattice after another period of a few hundred 
trap cycles, and then the lattice slowly damps over a further period of about 1000 
trap cycles. Lobo et al. made an approximate estimate of the total energy transferred 
irreversibly out of the condensate, and by assuming equipartition over the available 
modes, obtained a temperature of the thermal cloud which they assumed was respon- 
sible for the dissipation. More recently, Wright et al. [207] have treated this stirring 
problem using a TWPGPE approach. The initial vacuum noise gives an irreducible 
mechanism for seeding the dynamical instability, and their method ensures parti- 
cle number and rotating frame energy are conserved to very high accuracy over the 
length of the simulation. Furthermore the basis choice and numerical method they 
use is free from grid method and boundary artifacts such as aliasing and spurious 
damping at high momenta, so that any thermalisation and damping observed can 
be unambiguously attributed to the intrinsic field theory, rather than numerical arti- 
fact. The authors considered systems in 'pancake' traps, which are effectively two- 
dimensional. A typical final state, after the system has been subjected to a constantly 
rotating elliptical perturbation for 3000 trap cycles, is shown in figure [24)^a). 

This treatment allows a detailed and quantitative description of the thermalisation 
of the condensate. The thermal cloud created is initially located primarily in an outer 
annulus, and quickly obtains a classical moment of inertia, while the central region 
of the field is irrotational until penetrated by vortices. The temperature and chemical 
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potential of the thermal cloud are obtained by a self-consistent fit, and good agree- 
ment is obtained to an analytic estimate. The Penrose-Onsager criterion (see section 
3.2.4[ ) for identifying the condensate component fails in this system, due to the com- 



plex phase and amplitude structure associated with the vortex array. As an alternative 
method of characterising the coherence properties of the system, the authors examine 
the local behaviour of the temporal power spectra of the classical field about time to, 
namely 



H{x, oj; to) 



■ro+T/2 

iAc(x, t)e-''^'dt 

T/2 



(149) 



The sampling period T is chosen to be a number of trap cycles, and is long com- 
pared to the timescales characterising the phase evolution of the condensate (fi/fi), 
but short compared to the relaxation time of the field. Spectra such as shown in figure 



24 (b) are obtained. The spectrum from the central region of the field has a prominent 
and narrow peak at w w 10a»^, which is interpreted as the condensate eigenvalue. At 
larger radii (r > 9ro) the spectrum is broadened and is approximately that of the 
non-interacting gas. The local correlation times obtained from these data (by Fourier 
transform of H) allow an unambiguous distinction to be made between superfluid 
turbulence and thermal gas (which has a much shorter correlation time). We note 
that a feature of the 2D system is that the final state is not a regular Abrikosov lat- 
tice, but instead is a spatially disordered vortex liquid state. This can be interpreted 
as a thermally excited vortex lattice, and indeed both the simulations and an analytic 
prediction given in [207] show that a considerable amount of thermal energy is gen- 
erated in the stirring process, and that therefore the final condensate state must have 
considerable thermal excitation. The final state of the simulations is consistent with 
the condensate being in thermal and rotational equilibrium with the thermal cloud. 



4.7 Quantum statistical effects in Superchemistry 

The field of superchemistry was defined by Heinzen et al. [92] as "the coherent stim- 
ulation of chemical reactions via macroscopic occupation of a quantum state by a 
bosonic chemical species". Truncated Wigner simulations were used by Olsen et 
al. [151-155] to investigate superchemistry based on photo-association of trapped 
BECs into molecular dimers. The atom-molecule coupling occurs through a Raman 
two-photon transition for which the interaction Hamiltonian can be written as [152] 

Ant = ^ J d\x{^)(^\\^)^m'{^) - f,{x)^U^)) 

+i f t/^Xa(x)(«^^.(x)iAm(x)-«Am-(x)«A;„(x)) (150) 
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Figure 25. Comparison of atom number in photoassociation dynamics for Fock (solid line), coherent (dash-dot) 
and crescent quantum states of the initial atomic BEC. Reproduced from [153] © 2004 by The American Physical 

Society. 



where i^oCx) is the atomic field, ^„,'{x) the excited molecular field, and i^mCx) is the 
molecular ground state field, ;^f(x) is the Rabi frequency of the transition from atoms 
to excited molecules and Q(x) is the Rabi frequency for the transition from excited 
to ground molecular states. 

A notable feature of this work was the departure from the mean field predictions 
(see also [98]). The quadratic dependence of ( 150| ) on atomic fields combined with 
the relatively short timescale of the atom-molecule transition combine to make the 
process highly sensitive to quantum statistical effects such as squeezing (e.g. see 
figure[3]for comparison of some different quantum states for the condensate). Results 
of those studies demonstrated a regime where the quantum statistics of the atomic 
condensate play a crucial role in superchemistry dynamics and confirmed that photo- 
association of a trapped BEC into molecules provides a signature of the quantum 
state of the BEC. This can be seen in figure 25 where photoassociation dynamics are 
compared for different quantum states of the initial BEC. 



4.8 The quantum linewidth of an atom laser 

Johnsson et al. [112] used the TWA to model the process of weak outcoupling from 
a trapped Bose-Einstein condensate and to determine the linewidth of the outcoupled 
beam. Semi-classical analysis predicts the linewidth will be essentially Fourier lim- 
ited, scaling as the inverse of the outcoupling time [111]. Determining the quantum 
linewidth requires a multimode quantum theory of the atom laser outcoupler which 
the was implemented using the TWA. For Raman transition based outcouplers the 
primary source of linewidth broadening in the weak outcoupling regime comes from 
phase fluctuations of the condensate. The main source of phase fluctuations arises 
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from the nonlinear interactions in the condensate which convert number fluctuations 
into phase fluctuations [215]. A simple estimate for the quantum linewidth, A£, of 
the output coupled atom laser beam from condensate containing A'^o atoms in a har- 
monic trap is given by 

^E = ^AA^o, (151) 



where //tf is the Thomas-Fermi chemical potential (see equation ( 103 1). If the quan- 
tum state of the BEC is approximately a coherent state then the fluctuations in con- 
densate number are Poissonian, i.e. AA^o - ^JNq, and an analytic expression for A£ 
can be calculated. 

In [112] one and two dimensional truncated Wigner simulations were used to ob- 
tain the atom laser linewidth as a function of output coupling time, and these results 
were compared against GPE simulations. For short times the linewidth was found to 
be inversely proportional to the output coupling time, a feature adequately described 
by the GPE. However, on longer timescales the linewidth predictions of the two the- 
ories differed: the truncated Wigner simulations plateaued towards the Poissonian- 



limit ( 151 1, whereas the GPE continued to narrow. 



5 The stochastic projected Gross-Pitaevskii equation 

The stochastic projected Gross-Pitaevskii equation (SPGPE) is a truncated Wigner 
theory of Bose gases which takes into account the interactions between the atoms in 
the c-field region and the I region. The theory is valid for sufficiently large systems 
for temperatures from about O.STc to just above the BEC transition at T = Tc when 
the I region contains many weakly populated thermal modes. For a trapped system 
with largest trap frequency to the condition fiu «c hgT must be satisfied; in this sense 
it is a high temperature theory, extending the PGPE theory treated in section 3 and 
complimenting the low temperature TWPGPE approach described in section 4 An 
important point of difference is that the SPGPE theory is a grand canonical theory, 
in contrast to the microcanonical approaches described in previous sections. In its 
simplest implementation the theory is parameterized by the (in general time depen- 
dent) temperature T and chemical potential ji of the thermal reservoir comprised of 
the thermally occupied modes contained in the I region. Thus, in contrast to micro- 
canonical approaches for which temperature must be determined a posteriori, the 
SPGPE formalism allows direct control of the temperature of the interacting system. 
In general the I and C regions may be out of equilibrium (such as during quench 
cooling towards Tc). The dynamics of the condensation process are particularly well 
suited to treatment with this theory. 
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5.1 Formalism 

The treatment of thermal processes using stochastic methods has a long history, be- 
ginning with the theory of Brownian motion [66, 128]. The theory of open quantum 
systems couples the modes of interest in the quantum system to a reservoir [78]. The 
precise details of the derivation can be found in [31,32,73,74]. As a full treatment is 
rather lengthy, here we will briefly outline the development of the theory. 

The essential conditions for a treatment of the degenerate Bose gas with minimal 
complexity are i) the system and the reservoir are uncorrected, ii) the reservoir is 
in local equilibrium, described by the semi-classical Bose-Einstein distribution. In 
many systems of interest these conditions can be readily satisfied by appropriately 
choosing the cutoff" energy ecut. In essence, the stochastic PGPE theory extends the 
PGPE theory by including the coupling of C region atoms with atoms in I which form 
a grand-canonical reservoir. The coupling generates additional damping and stochas- 
tic terms in the PGPE, which necessarily takes the form of a stochastic differential 
equation. 

The assumption of local equilibrium for the I region is convenient, but does not 
pose a fundamental limitation of the formalism. In principle it is possible to derive a 
quantum Boltzmann-Iike kinetic equation for the particles in I coupled to a stochastic 
c-field equation for the particle in C. However, this would result in further compu- 
tational complexity which is not necessary for a broad range of applications. Such a 
description remains a goal for the future, and would result in a near complete model 
of Bose gas dynamics at high temperature. 

5.7.7 Background 

The theory of finite temperature BEC dynamics of Zaremba, Nikuni, and Grif- 
fin [219], which was developed along the hnes of the two fluid theory of superfluid 
helium, provided the foundations of generalized GPE theory from a hydrodynamic 
point of view. The essential assumption of the theory is that atoms enter and leave the 
condensate so as to enforce local energy and momentum conservation. The resulting 
description takes the form of a finite temperature GPE for the condensate, coupled 
to a Boltzmann equation for the non-condensate. The theory has the advantage that 
the condensate and non-condensate are described on an equal footing, making the 
description of coupled dynamics of thermal cloud and condensate tractable. While 
being intuitively appealing and a providing a powerful approach for extending zero 
temperature mean-field theory, enforcing local conservation of energy and momen- 
tum also has disadvantages which are removed by the SPGPE approach. Firstly, as 
a mean field theory it cannot be valid near the BEC phase transition where the con- 
densate is relatively small or non-existent. Furthermore, a locally conserving theory 
neglects the non-locality of quantum mechanics, which plays a fundamental role in 
determining the dynamics of atoms stimulated into a highly occupied field. 
The SGPE formulation of Stoof is closely related to the SPGPE theory presented 
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here: the formulation leads to a stochastic differential equation for the condensate 
which is driven by a noise term associated with condensate growth. However, there 
are two important differences. Firstly, in [201] the reservoir is chosen to contain 
all modes with energy in excess of the chemical potential ji. The stochastic GP- 
equation so obtained involves self-energy functions, necessitating a many-body T- 
matrix treatment of interactions. As well as being difficult to implement numerically, 
the c-field of the theory only describes the condensate and few very low energy ex- 
citations. As discussed in section 3.1.2 in the vicinity of the transition there are 
typically 10-^ - lO'* degenerate modes warranting a c-field treatment. Secondly, the 
approach neglects scattering terms which conserve population but transfer energy 
from the reservoir to the c-field region. The inclusion of these terms in the SPGPE 
stems from the explicit use of a high energy projector, which renders them finite and 
tractable. 



Stochastic PGPE. The Quantum Kinetic theory of Bose-Einstein condensation [49, 
76, 77] has been shown to provide a good description of the process of conden- 
sate formation [75,79, 123] at the level of condensate population dynamics provided 
condensation occurs into the absolute ground state of the system. Recent work by 
Gardiner et al. [30-32,73,74] developed the SPGPE theory to explicitly include a 
high energy cutoff, thus unifying the PGPE theory with the reservoir theory of high 
temperature BECs. 

Since its inception, the theory of finite temperature BECs has presented many chal- 
lenges and a rich body of work has accumulated [176]. Several key problems are: 
the consistent treatment of two body interactions [101, 143], the description of ther- 
mal cloud coupling and dynamics [201,219], and the unification of GPE approaches 
with dissipative finite temperature phenomenology [40,209]. In these respects, the 
PSGPE theory has some notable computational and physical advantages which we 
briefly summarize: 

(i) Consistent UV-cutoff: The imposition of a high energy cutoff using the methods 
of section [3] imposes a consistent cutoff, even for trapped systems. As noted in 
section [T. 1.1 1 if the cutoff is only imposed in momentum space the precise wave- 
length for the cutoff is position dependent. The PGPE formalism addresses this 
problem by imposing a cutoff in the single particle basis which approximately 
diagonalises the many-body Hamiltonian at high energies. 

(ii) Two body T-matrix: By imposing the cutoff at high energies the need to use the 
many body J-matrix to describe scattering is eliminated, requiring only the two 
body T-matrix description of S-wave scattering: T{0) AnTp-ajm. 

(iii) Nonlocal description of condensate growth — Fundamentally the theory is non- 
local: atoms which leave the high energy cloud enter the c-field region so as to 
minimize the difference between the chemical potential of the high energy region 
and the GP operator acting on the c-field. Beyond hydrodynamic effects are thus 
explicitly included. 



78 

(iv) Scattering terms: Imposing a cutoff at high energies allows the so-called scatter- 
ing terms — reservoir interactions that do not directly change the populations of 
the reservoir or c-field — to be consistently included in the theory. Analogous 
terms proved to play an important role in the Quantum Kinetic theory description 
of condensate growth [130]. In the SPGPE theory the scattering terms couple to 
dynamical excitations in the c-field region. 

(v) Valid at the BEC transition: A notable feature of the SPGPE is that it is particu- 
larly well suited for dynamical studies in the regime T ~ Tc since the conditions 
of validity are high temperature (hoj <sc keT) and moderate occupation of modes, 
both of which are readily satisfied near the BEC transition. 

(vi) Reservoir dynamics: In principle the dynamics of the incoherent region can be 
treated with a Quantum Boltzmann equation with little additional formalism pro- 
vided the region remains in approximate local equilibrium. 

(vii) Consistent mean field treatment of dissipation The mean field theory recovered 
by setting all noises to zero gives a GP equation of motion with extra dissipative 
terms. The dissipation evolves the c-field to a ground state with the chemical 
potential of the reservoir. In contrast to phenomenological approaches [40, 209] , 
the dissipation rate arising in the equation of motion is the physical rate of the 
theory which can be derived from a Boltzmann integraQ[32]. 

5.7.2 The system and its separation 

We again consider a dilute Bose gas held in a trapping potential defined by ([3]), but 
extend the PGPE description of section [3] to consider the coupling of the c-field re- 
gion to the I region. To accommodate applications such as the formation of vortex 



lattices in BECs (discussed below in section 5.4 1, we include the possibility that the 
I region may be rotating. This requires either that the trapping potential is axially 
symmetric, or time-independent in a rotating frame of reference (corresponding to 
elliptical stirring at a constant angular frequency). For simplicity we restrict our at- 
tention to systems where either i) the I region is stationary in the laboratory frame, 
or ii) the I region is rotating in an axially symmetric trap. For the latter case the the- 
ory is conveniently formulated in the frame rotating at the frequency of the I region, 
which we denote by Q. Choosing the symmetry axis of the system to be the z-axis, 
the single particle Hamiltonian for the system transformed to the rotating frame is 



^sp = Ho = + Voix) - aL„ (152) 

2m 



where = -ih(xdy - yd^) is the z-component of the angular momentum operator. 
It has been shown that when the incoherent region is in rotational equilibrium in a 



The difference between ttie effective dissipation rates in the two descriptions is usually small for weak dissipation, 
but the equations of motion are formally distinct. See reference [31] for a more detailed comparison of dissipative 
mean field approaches to vortex formation in BECs. 
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harmonic trap the theory is modified by transforming the c-field description to the 
rotating frame and making the replacement ojr ^ a)_i_, where 



a)±_ 



^0)1 - 02, (153) 



in the dissipation rates of the SPGPE [32]. We can thus treat the rotating case in 
the formalism by including the effects of rotating frame transformation in the c- 
field description (152]), parameterized by Q. We return to rotating systems in more 



detail in section 5.4.2[ but in what follows the formalism applies to systems that are 



either non-rotating and in general non-axisymmetric, or rotating and axi-symmetric 

(CU^. ^ = (jJr). 

As described in section |2.2[ the field operator for the full system is decomposed 
into a c-field and an incoherent field. The SGPE takes the form of an equation of 
motion for the c-field i/'cCx) with terms arising from interaction with the incoherent 
region I. The two regions are treated using different approximations. 

5.7. i Treatment of the incoherent region 

The local equilibrium assumption for the state of the incoherent region allows all 
higher order correlation functions arising in the theory to be factorized into products 
of second order correlation functions. At this level of approximation the essential 
reservoir interaction physics can be reduced to functions of the single particle Wigner 
function for the incoherent region 

Fi(x, K) = J d\' ((A!(x + x72)«Ai(x - x72)) ^>'''"', (154) 



previously introduced in section 3.2.6 Introducing the variables u = (x -i- x')/2, 



V = x' - X, we can write 



((A!(x')«^i(x, t)> = (.a!(u + v/2)«Ai(u - v/2, t)), 

« f J^Kfi(u,K)e-'K™("''^)\ (155) 

<^i(x')<a!(x,t)>« r j3K[l+Fi(u,K)]e'*-^'«)r^ (156) 
where phase space integration over the incoherent region Q.\ constrains the coordi- 



nates to satisfy ?iaj(x, K) > ^cut (see 124 and 127 1 and the energy in the frame rotating 



with angular frequency vector Q = Qz has the semi-classical form 

Hojix, K) = — nil-{xxK) + Vo{x). (157) 

2m 
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The dissipation rates of the theory are time-integrated products of such functions, 
as can be seen from ( |171| ) and ( |174| l. 

Semi-classical Bose-Einstein distribution. For many applications the I region may be 
described by a semi-classical Bose-Einstein distribution 

" uf. i ^/^. n i ' ^^^^^ 

exp [(nw(x, K) - ^)lkBT] - 1 

where we note that for high energies where this is assumed to apply, the gas is always 
rapidly thermalized and interactions with the c-field region are a small correction to 
the single particle energy ( 157 1. This form has been used to evaluate the dissipation 
rates of the SPGPE theory [32]. 

The properties of the ideal gas description of ( 158| ) and ( 157 1, including the effect 
of the cutoff, can be expressed in terms of the incomplete Bose-Einstein function 
defined as [32] 



1 /^oo 

r(y) J, 



(159) 

ZZ^ r(v, yl) 
1=1 r(y) ' 

where T{v,x) = dy is the incomplete Gamma function. In analogy with 

the reduction to an ordinary Gamma function r(v, 0) - r(v), we have gy{z,G) = 
9viz) = HJ^i z' /l^, reducing to the ordinary Bose-Einstein function. We then find for 
the I region density, for example 



"iW = I 7;rTJ^i(x, k) 



where ffKcutix)^ /2m = max{ecut - V'j_(x),0) and we have introduced the effective 
potential 



m 



W(x) - --(ojy + (161) 

which accounts for the transformation to the rotating frame ([153 1. Setting e^ut = 0, we 
recover the standard form for the semi-classical particle density of the ideal gas [43] . 
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The total number in I region is given by 



M - 03 (e^^ySecut) /OS?^a))^ (162) 



where cb = {(jo^aPi)^^^ is the geometric mean frequency in the rotating fame. In this 
way the usual semi-classical expressions can be generalized to include a cutoff in 
terms of the incomplete Bose-Einstein function. 



5.7.4 Treatment of the c-field region: deriving the equation of motion 

The standard procedure of phase space methods for open systems involves deriv- 
ing a master equation for the reduced system by eliminating the reservoir degrees of 
freedom. The master equation may then be mapped to a generalized Fokker-Planck 
equation (FPE) of motion for a quasi-probability distribution, such as the Wigner 



representation, by making use of operator correspondences (e.g. (52)-(55l). Pro- 
vided the FPE contains derivatives which are at most second order (representing 
diffusion), an equivalent stochastic differential equation may be found which can be 
conveniently simulated numerically. The truncated Wigner approximation involves 



neglecting third order terms in the FPE (see (58 1), and in the context of the SPGPE 
theory the truncated Wigner approximation reduces the FPE to second order, allow- 
ing a formulation of the problem in terms of stochastic differential equations. 



Validity of the SPGPE treatment of the c-field. In addition to treating the I region semi- 
classically, the SPGPE formalism makes the truncated Wigner approximation which 
requires that the modes under consideration are highly occupied. An approximate 
master equation that can be mapped to a stochastic differential equation is then 
obtained by truncating the interaction between C and I at first order in powers of 
hto/ksT, where oj = max{w,} is the largest oscillator frequency of the system. 

A feature of the the high temperature theory is that the dissipation arising from 
the reservoir coupling acts to smooth out sharp phase space structure that would 
otherwise generate significant third order term corrections [223]. In this sense the 
high temperature Bose gas is particularly well suited to treatment using the truncated 
Wigner method. 
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5.1.5 Stochastic projected Gross-Pitaevskii equation 

In the rotating frame the full non-local form of the SPGPE is given by the following 
stochastic differential equation in Stratonovich form 

(5)#c(x, t) = !Pc| - ^Lc«Ac(x)c?f (163a) 
+^(M-Lc)M^)dt + dWG(x,t) (163b) 



+ f d'x' M(x - x') ^^^—^^^-^iPc{x)dt + /tAc(x)t/WM(x, \. (163c) 



The first line of the SPGPE (163ai describes Hamiltonian evolution according to 



the POPE introduced in section 2.3.6 and developed in section [3] The projector 
appears as a natural consequence of formally imposing a high energy cutoff in the 
definition of the c-field region. The operator Lq is the Hamiltonian evolution operator 



for the c-field region defined via equations (65 1 and ([67]). Its explicit form is 



Lc<Ac(x) = (//sp + m|(Ac(x, tf) (^c(x). (164) 

The second line of the SPGPE ( |163b[ ) is directly responsible for condensate growth 
from scattering between two I region atoms as illustrated in figure 26 a). The yu and 
T that arise are respectively the chemical potential and temperature of the thermal 
reservoir of particles in the I region. The quantity G(x) is a spatially dependent colli- 
sion rate, specified by a quantum Boltzmann integral over the I region as discussed in 



more detail in section section 5.2.1 below. The complex noise associated with growth 
is dWcix' , t) and satisfies 

(JWg(x, t)dWG{x:,t)) = 2G(x)dc(x,x')dt, (165) 
{dWcix, t)dWG{x',t)) = {dWc(x, t)dWc{x', t)) = 0. (166) 



The third line of the SPGPE ( 163c| ) represents number conserving scattering pro- 



cesses between atoms in C and I and provides a mechanism for energy transfer be- 
tween the two regions as illustrated in figure[26^b). This couples to the divergence of 
the c-field region current given by 

7c(x) ^ ^([V.Ac(x)]'Ac(x) - (Ac(x)V«Ac(x)) - (Q X X) |(Ac(x)P, (167) 

where the last line is a rigid-body rotation term arising from the transformation to 
the rotating frame. The limit 7c(x) = gives the laboratory frame velocity field 
V = Q X X, corresponding to the irrotational system mimicking rigid body rotation. 
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(a) Growth 




Figure 26. Schematic of the processes arising from the interactions between the c-field and incoherent regions. In 
(a) two c-tield region atoms collide, with a significant fraction of the collision energy transfeired to one of the 
atoms, with the other atom passing into the c-field region. The time-reversed process also occurs. In (b) a c-field 
region atom collides with a incoherent region atom with no change in c-field region population. 



The rate function M(x - x') is specified by a second quantum Boltzmann integral, 
and is discussed in detail below in section section [F.2.21 
The real noise (/Wm(x, t) associated with scattering is specified by 



(dWMix, mWui^a', t)) - 2M (x - x') dt. 



(168) 



Grand canonical equilibrium. As described in section 3.2.4[ the PGPE provides a 
means to sample the microcanonical ensemble of equilibrium states. By including 
interactions with the incoherent region we have arrived at a grand canonical descrip- 
tion, parameterized by the chemical potential and temperature of I. Irrespective of 
the form of G(x) and M(x), the SPGPE evolves the system to the grand canonical 
equiUbrium distribution 



Ws exp 



uNc - He 
kuT 



(169) 



corresponding to the density matrix 



Ps oc exp 



liNc - He 
kuT 



(170) 



in the truncated Wigner approximation. Once the c-field reaches equilibrium single 
trajectories may be used to sample the grand canonical ensemble. 
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Figure 27. Dependence of the growth rate y jl72) on {i for the choice Eciit = 3/J. The full rate jl72| (solid line) is 
compared with the logarithmic term (dashed line). The two points are calculated for N = 10* ^ Rb atoms in a trap 
with geometric mean frequency = 2ff X 25Hz using the ideal gas relation for T,.. The chemical potentials are 
estimated from fi a 3Pioj/2 aXT = T^, and fj » jUtf(^o) at T = 0.57^. 



5.2 Growth and scattering in the SPGPE 



We now discuss the properties of the dissipative terms in the SPGPE ( 163 1. We give 
the explicit form of the rate functions G(x) and M(x), the regimes under which they 
may be evaluated in closed form, and discuss details of their physical interpretation. 

5.2.1 Growth terns 

Growth rate. The explicit form of the growth rate is [32, 74] 



G(x) = — — I I I d'Kid^K2d^K3F{x,Ki)Fix,K2) 



{Inrn- jjjsii 

X[l+F(x,K3)]Ai23(0,0), (171) 

where Ai23(k, e) = d{K\ + K2 - K3 - k)6(coi + aj2 - - e/h) conserves energy 
and momentum during the collision. The rate G(x) can be calculated in the regime 
where the I region is quasi-static, and well-approximated by an ideal semi-classical 
Bose-Einstein distribution ( 158 1 with slowly varying and T{t). For the inner 
spatial region of a harmonic trap satisfying V{x) < IEr/S we find that G(x) = y is 
independent of position with 

y = y„| [in (1 - + ^.2/^(/^--) J] ^'■/^a'-2^e„,) (cD[/(^-«\ 1, r + l]f |, (172) 



where yo = 4m(akBT)^ /nH^ and (^[x,y,z] is the Lerch transcendent. Outside this 
region there is a weak spatial dependence which can be neglected for most purposes 
[32]. The bare rate, 7o> has been used as an estimate in the literature, often requiring a 
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'fudge factor' (usually chosen as ~ 3) to obtain a rate that gives physically reasonable 



damping times. In figure 27 we show y for a fixed choice of tcut = 3//. The full picture 



is more complicated than shown in figure 27 because the choice fcut = 3/^ would not 
be appropriate near Jc, but instead a much higher fcut would be necessary. In practice, 
the more accurate form typically increases the ratio y/yo by a factor which is in the 
range 1-10. 

Dissipative dynamics of condensate growth. By neglecting the scattering and noise terms 
in the SPGPE ( |163| ) it is possible to show that 



d{Hc-fiNc) 
dt 



4J 



d\ |(//-Lc)<Ac(x,Or, 



(173) 



where we used the approximation G(x) ~ y d& given by equation ( 172 1. As the RHS 



of equation ( 173 1 is a strictly non-positive term, we can see that the growth term 
acts to minimize the effective grand-canonical Hamiltonian = tic - l^^c- The 



equilibrium solution is the ground state of the PGPE (68 1 with chemical potential /u 



The growth terms describe the Bose-stimulated transfer of particles between the C 
and I regions during two body collisions. 

5.2.2 Scattering terms 



Scattering rate. The rate function for the scattering term of the SPGPE ( 163 1 is most 



easily calculated by transforming to momentum space M(x) - J<i'^k 
where we find 



M(k), 



M(k) = -^ rr J^KiJ^Kj Ai2(k,0)f(u,Ki)[l+f(u,K2)], (174) 

with Ai2(k, e) = 6{a>i + a)2 - e/fi)d{Ki -i- K2 - k). It has been shown that to a good 
approximation this expression is independent of u [32,74], and for this reason we 
suppress the argument in the present definitions. 

Using the same approximation of a quasi-static thermal cloud as used in the calcu- 
lation of the rate G(x) in the previous section we find that 



M(k) 



\6na ksT e> 



M 



{2nffi\\i\ (ef^(E!<-^l) _ if (27r)3|kr 



(175) 



and thus 



M(x) 



(2)r)' J 



-/kx 



|k| 



(176) 



so that M(x) is a spatially dependent function over the whole C region. At first glance 
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Figure 28. Scattering potential for tlie breatliing mode Sn i (r, t) of a Tliomas-Fermi condensate, (a) Tlie radial 
sliape of tlie mode f\(rlRif) (see text), (b) The scattering potential for a breathing mode excitation with amplitude 
A = 0.1 (solid line, shown at ; = 0) and the harmonic trapping potential (dashed line). We have used Ecut = 3jU and 
// » fiTp(No) for T = OJTc and have estimated Nq using the ideal gas relation for 10* ^'Rb atoms in a trap with 
radial frequency oj^ = 2nx 25Hz. For these parameters VM.iir) is small compared to the harmonic trap for 

significant radii (shown inset). 



this term appears somewhat pathological, but well defined results are obtained since 
M(x) is convolved with functions of condensate band fields. Such functions are both 



UV- and IR- cutoff, giving a finite result for the convolution in ( 163c I 



Effect of scattering on hydrodynamic collective modes. To gain some physical insight 
into the nature of the scattering we note that the evolution according to the determin- 



istic part of ( 163c i can be written as a real effective potential 



ifi 



5«Ac(x,0 



dt 



-:Pc{Vm(x,0'Ac(x,0}: 



(177) 



M 



where 



VmCx, - - J d\' M (x - x') ^ V • 7c(x'). (178) 



To first approximation we can neglect all dissipation as relatively weak corrections 



to the PGPE evolution ( 163a i, giving the continuity equation V • 7c(x) ~ -dncix)/dt 
for the c-field region and 



Vm(x, t) 



M(x-x') 



dnc(x') 
kTf dt 



(179) 



Thus the scattering term generates an effective potential from dynamical density 
fluctuations in the c-field region. As a specific example we consider excitations of 
a system with spherical symmetry, since the momentum distribution will have the 
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same symmetry as the scattering kernel ( 176| ). We can then evaluate the effective 
potential for collective modes in the hydrodynamic and Thomas-Fermi approxima- 
tions. For a spherically symmetric trap the condensate wavefunction in the Thomas- 
Fermi approximation is «tf(x) = (jj.tf/u)(1 - (r/7?TF)^), with Thomas-Fermi radius 
Rjp = yjlidjp/mcoj. The density profile for spherically symmetric modes [204] with 
amplitude A can be written as 



Auyf 

6n„(x,t) = sm(aj„t)f„(r/Rjf), 



(180) 



where the radial form is 

(„.i/yV(o,i/2)(2,2 

/«(0) = ' ' 



given by the 

l)6'(l - x) and 9{x) 
1 the peak density of the excitation is AyUxp/M- The modes have fre- 



Jacobi polynomials fn{x) = 
is the unit step function. Since 



quencies a>n = cOr ^lln^ + 3n, f or ex ampl e the breathing mode has frequency 



cl>i = VSctif. Evaluating equations ( 179 1 and ( 176 1 leads to 



(4aRl 



yM,„(x, t) = Anaj„ cos ioj„t) f j ^-^^_-^_^F„(r/7?TF), (181) 



where a,- = ^|W|mu^r is the radial harmonic oscillator length, and 



sin(itx) r' 

F„{x)= dk— dy ysm(ky)f„{y). (182) 

Jo kx Jo 



In figure 28 we show the effect of scattering on the breathing mode. At t = the 



density modulation ( 180l vanishes, but has maximal rate of change, reflected by the 
phase of (TSTJ relative to ( 180l. The radial shape of the mode fiir/Rjp) is shown in 
figure 28 1, corresponding to outward flow at f = 0. The potential, shown in figure [28|3 
(at t = 0), generates damping of the excitation by imposing an additional potential 
gradient that acts to oppose the outward flow. The same qualitative result holds for 
higher modes, where VM,nir) is found to have the same overall shape as fnir/rjp) and 
a relative phase so as to oppose the excitations with an additional potential gradient. 
From figure 28 d it is clear that VM,i(r) can be a significant correction to the bare 



trapping potential near the center of the trap, but is unimportant near 7?tf (see inset). 



5.3 Simple growth SPGPE 



The scattering term of the SPGPE ( 163c i does not alter the population of the conden- 
sate band directly, nor does it affect the grand canonical equilibrium of the c-field. It 
is generally expected to be less important to the c-field dynamics in comparison to 
the growth term ( 163b[ ), which directly describe the dominant collision processes re- 
sulting in Bose-Einstein condensation. Combined with the difficulty in its numerical 
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Figure 29. Spontaneous formation of vortices during Bose-Einstein condensation, a, 200-//m-square expansion 
images of BECs created in a liarmonic trap, sliowing single vortices (left, centre) and two vortices (right), b, c, 
Sample simulation results from evaporative cooling in a harmonic trap, showing in-trap integrated coluimi densities 
along z (in b) and associated phase profiles in the z = plane (in c), with vortices indicated by crosses and circles at 
±2n phase windings, d. Left image: 70-yum-square phase-contrast experimental image of a EEC in a toroidal trap. 
Remaining images: vortices in 200-//m-square expansion images of BECs created in the toroidal trap, e, f. 
Simulations of BEC growth in the toroidal trap show vortices (as in b,c) and persistent currents. Reproduced 

from [212] 

implementation arising from the nonlocality of the deterministic part and the multi- 
pUcative nature of the associated noise, it seems reasonable to neglect it in the first 
instance. This results in the simple growth SPGPE 



#c(x,0 ^Pc\-Ucil^c('!iJ)dt + -P^Qi - Lc)il^ci'ii,t)dt + dWy{x,t) \ , (183) 
n kbT 



where 



{dW*Jx,t)dWy{x',t)) - 2ydc{x,x')dt. 



(184) 



The numerical implementation of this simplified form is relatively straightforward, 
being only somewhat more complicated than the PGPE. As the only noise term is 
additive, the simple growth SPGPE can be integrated using high-order algorithms, 
such as a modified fourth-order Runge-Kutta algorithm [141]. 



5.4 Applications to the dynamics of partially condensed Bose gases 

Background. The first stochastic Gross-Pitaevskii treatment of Bose gases was devel- 
oped by Stoof [201] using a functional integral formulation of the Keldysh method. 
The first application of SGPE theory was carried out by Stoof and Bijlsma [203] who 
used the SGPE theory developed in [201] to study finite temperature dynamics of a 



89 



one dimensional Bose gas. Focussing on the scenario of growth into a one dimen- 
sional dimple trap the authors studied reversible condensate formation (see [196]), 
and the frequencies and damping rates of collective modes. The SGPE theory has 
also been used in conjunction with variational techniques to study finite temperature 
collective excitations [64,203], and dissipative vortex dynamics [63]. More recently, 
Proukakis et al. [174] have used the method to investigate quasicondensate growth 
into a one dimensional dimple trap. For a deep dimple the dynamics were found to 
involve shock-wave propagation in addition to quasicondensate formation. 

Applications of the SPGPE theory are still relatively few in number. The cur- 
rent authors have only recently completed the implementation of the simple growth 
SPGPE for a harmonic trap in two dimensions with rotation, and in three dimensions 
without. However, even the dissipative mean field equation obtained by neglecting 
all noise terms in the SPGPE has given significant insight into the dynamics of fi- 
nite temperature BECs. Here we briefly review the applications that have appeared 
to date. 



5.4.1 Spontaneous vortex formation during Bose-Einstein condensation 

The formation of topological defects in symmetry-breaking phase transitions has 
been a topic of interest in both cosmological [120] and condensed matter [222] 
scenarios. Until recently quantitative models of the formation dynamics of trapped 
Bose-Einstein condensates have only been studied at the level of populations with- 
out allowing for the possibility of topological excitations [20,48,49,75,79, 123]. The 
possibility of the spontaneous formation of vortices in the growth of a trapped Bose- 
Einstein condensate has been suggested previously by Anglin and Zurek [6] and 
Svistunov [206], but until recently had not been confirmed by experiment and the 
possibility has been the subject of much speculation [202]. Recent work by Weiler et 
al. [212] has reported the observation of spontaneous vortex formation in the growth 
of a trapped BEC, and compared the statistics of formation with predictions of the 
simple growth SPGPE ([T83]). 



Experiment. The experiments of Weiler et al. [212] evaporatively cooled a dilute 
gas of ^^Rb atoms from slightly above Tc in both an oblate harmonic trap and a 
toroidal trap formed by the addition of a Gaussian barrier from a tightly-focussed 
blue-detuned laser beam along the symmetry axis. The oblate nature of the trapping 
potential resulted in an energy penalty for vortices not aligned with the symmetry 
axis, hence improving the fidelity of vortex detection. After condensate formation 
and sufficiently long time-of-fiight expansion, vortex cores were observed with a 
probability ranging from 20-60%. 



Theory. The results of Weiler et al. [212] were modelled using the simple growth 
SPGPE (183 1 matched to the growth of the condensate number in the experiments 
by ramping the chemical potential and temperature of the thermal cloud. For each 
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Time= 0.13 s 



N„ 



352 



Time = 0.45 s 
26205 




Time = 0.67 s 
N^= 113311 





Time = 1 .57 s 
348327 



Figure 30. BEC growth dynamics, a-d. Four snapshots during the simulated growth of a BEC showing isodensity 
surfaces (in light red) in a three-dimensional rendering. Vortex cores of opposite charges about the z axis are 
indicated as magenta and cyan lines. The corresponding times are a, 0.13 s; b, 0.45 s; c, 0.67 s; d, 1.57 s, where 
/ =0 s is the time when the quench is initiated in the simulation. Reproduced from [212] 

experimental evaporative cooling ramp an ensemble of ~ 300 SPGPE trajectories 
was computecj^ and the vortex observation statistics were compared with the exper- 
imental measurements. A comparison between experimental absorption images and 
column densities of representative SPGPE trajectories is shown in figure |29j with 
clear qualitative similarity. A comparison of the vortex observation statistics yielded 
quantitative agreement between experiment and the SPGPE theory. In figure 30 a 



time series of density iso-surfaces is shown for a particular trajectory. The emer- 
gence of the final BEC is seen to be a turbulent process, in this case resulting in a 
metastable vortex. 



5.4.2 Rotating Bose-Einstein condensation 

Bradley et al. [32] have used the simple growth equation ( 1 83 1 to model the dynamics 
of the formation of a rotating Bose-Einstein condensate [90] where it is possible for 
stable vortices to form during condensation. The following question arises: does a 



'See online supplementary information for [212] for a subset of trajectories showing possible outcomes for each 
scenario of the experiment. 
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vortex-free condensate form before it is penetrated by vortices, or does condensation 
proceed into a state with vortices already present? 

It is helpful to consider the single particle energy spectrum (in the rotating frame) 
of the cylindrically symmetric harmonic trap in three dimensions 

enhn ^ ficOr(2n + \l\ + 1) - JiQ.1 + haj^im + 1/2), (185) 

where n, 1, m are the radial, angular and axial quantum numbers. In the absence of in- 
teractions BEC will form in the ground state mode with energy eooo = Tito,- + fioD^Jl. 
Vortices arise from occupation of states with nonzero angular momentum and the 
positive angular momentum part of the spectrum behaves as fi{a)r - 0)l leading to 
near-degeneracy of positive angular momentum modes as O ^ ojr- States with angu- 
lar momentum (L,) = ?i/ > are increasingly easy to populate as the rate of rotation 
increases. 

For a rapidly rotating BEC transition, vortices can play a dominant role at all stages 
in the BEC growth process. It is then possible that atoms condense into a vortex-filled 
but spatially disordered state — a vortex liquid. 

SPGPE treatment of rotation. In our development of the SPGPE we have included the 
possibility that the thermal cloud occupying I may be in rotational equilibrium in a 
symmetric trap. The transformation has consequences for the choice of fcut and the 
basis of single particle states. In harmonic oscillator units the single particle modes 



corresponding to ( |185[ ) are also eigenstates of L^: 

<f>nim{r, e, z) = N„i^e''r^'^e-'"'^L„{r^)e-'"'^H^{z) (186) 

which is the appropriate basis for introducing a projection operator to separate the 
C and I regions. Imposing the cutoff in the rotating frame leads to a bias (in the 
direction of the rotation) for the angular momentum for the c-field region modes. In 
the limit of very fast rotation (Q. a)r), z-excitations are strongly suppressed and 
only the lowest Landau level (LLL) is contained in the c-field if tcut < ^<^z- 



In figure 31 we show a representative SPGPE simulation of the BEC transition 
in a rapidly rotating system (not in the LLL regime). An initial thermal state with 
no condensate present is evolved subject to a sudden change in the thermal cloud 
temperature and chemical potential consistent with a quench below the critical point. 
When the final temperature after the quench is comparatively low {T <«c Tc shown 
here) the dynamics show a phase of rapid growth into a vortex liquid, followed by 
a much slower ordering phase where the vortices assemble into a regular lattice. At 
higher temperatures {T ~ Tc), condensation into a vortex liquid still occurs but the 
ordering phase is frustrated by thermal fluctuations. 
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(a) uj,.t = 0.0 (b) Urt = 62.5 (c) uj^t = 125.3 
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Figure 3 1 . Rotating Bose-Einstein condensation: c-field region density for a single trajectory of the SPGPE (in the 
frame of the thermal cloud rotating at fio = 0.979(1),). (a) Initial state for 1.3 X 10'' ^'■'Rb atoms at Tq = 12nK, with 
fio = 0.5eooo. At time f = the non-condensate band is quenched toT = InK, fi = 3.5eooo, and H = fio preserving 
the rotation rate, (b)-(d) The c-field region undergoes rapid growth into a vortex liquid state, (e)-(f) At this low 
temperature the vortices then assemble into a regular Abrikosov lattice. Reproduced from [32] © 2008 by The 

American Physical Society. 

6 Conclusion 

In this review we have outUned a unified c-field theory for studying the dynamics 
and statistical mechanics of ultra-cold Bose gases. The various c-field approaches 
considered apply to situations as diverse as the zero-temperature quantum dynamics 
of colliding Bose-Eintein condensates, through to the effects of critical fluctuations 
at the condensation transition. Perhaps one of the most surprising aspect of these 
approaches is that underlying them is the well-known, and computationally obliging, 
Gross-Pitaevskii equation. Two essential adaptions to the Gross-Pitaevskii theory 
bring this power to describe quantum and thermal effects: stochastisation of the field 
and evolution equation, and projection onto the c-field region. 

Recently there has been increasing use of classical field techniques in the ultra- 
cold atom community, and particularly in the truncated Wigner approach to describe 
beyond-mean-field dynamics in BECs at zero temperature. Most of the classical field 
calculations performed at finite temperature have not made controlled use of projec- 
tors, other than the accidental projection intrinsic to the numerical representation. 
One of the consequences is that these results can often not be easily related to ex- 
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periment. In this review we have extolled the virtues of a consistent energy projector 
used to define the c-field region, and have shown that this allows a theory that can be 
applied to quantitatively describe experiment. 

The future development of c-field techniques is closely related to the direction in 
which many ultra-cold gas experiments are heading. We close this review by pointing 
to a few such areas under investigation: 

Full implementation ofSPGPE. To date the SPGPE scattering term has yet to be im- 
plemented numerically, primarily due to technical challenges. This term is expected 
to have an effect on strongly non-equilibrium scenarios such as condensate growth 
and unstable vortex dynamics, although it is currently not clear in what regimes it 
would give rise to measurable differences. Developing a formulation that can handle 
the dynamics of the incoherent region remains an important future challenge which 
would complete the theoretical description for most practical purposes. 

General interactions. Long-range interactions, such as dipole interactions, and strong 
s-wave interactions using Feshbach resonances are now routinely available in ex- 
periments. For these systems finite temperature and incoherent processes appear to 
become more important than in the weakly interacting s-wave case, and should be 
well-suited to a c-field description. 

Fermionic systems. The quantum statistics of bosons allows for the modes of the 
atomic field to be highly occupied, providing the coherence central to the c-field 
description of the ultra-cold Bose gas. In contrast, quantum statistics prohibits multi- 
ple fermions from occupying a single mode, and thus invalidating a c-field approach 
at face value. A fundamental question to address is whether there is some other path- 
way to providing a useful c-field description of Fermi gases. 
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Appendix A: Numerical technique for the harmonically trapped system 

In this appendix we overview a numerical method that allows an efficient and ac- 
curate solution of the projected Gross-Pitaevskii equation in a harmonic potential. 
In Appendix |B] we briefly discuss a method for the uniform gas and refer to refer- 
ences [23,27,32,207] for a more detailed discussion of implementation details and 
applications to rotating systems. 

A.l Numerical requirements 

The modes of the system are of central importance in the assumptions used to derive 
the various c-field methods presented in this review, and care must be taken in nu- 
merical implementations to ensure the modes are faithfully represented. Any useful 
simulation technique must satisfy the following requirements. 

(i) The space spanned by the modes of the simulation should match the c-field region 
as closely as possible. 

(ii) All modes in the c-field regime must be propagated accurately. 

The case we examine here is the PGPE equation for the harmonically trapped 
system, i.e. where 

VoCx) = ]^mo? [x^ + y^ + z^) . (Al) 

To simplify the discussion we have taken the harmonic trapping potential to be 
isotropic!^ and will not consider the perturbation potential 6V. We take the c-field 
region to be defined by an energy cutoff in the single particle basis, i.e. eigenstates 
of Ho = p^/2m + Vo(x) with energy less than ecut- 

A.2 Spectral representation of the PGPE 

For convenience, we write the PGPE in dimensionless units to simplify the discus- 
sion, and explicitly indicate all dimensionless quantities in this section by use of 
tildes. We do this by introducing a unit of distance xq = yjWJmo) and time fo = 1 1^- 
These choices immediately imply computational units for energy Eq = Hco and mo- 
mentum po - ^Jtima). So, for example, our dimensionless distance variable is defined 
as X = x/xo, dimensionless time, i = t/to, and c-field, i/'c = *Pc^^J^- The coefficient 
of the nonlinear term in the Gross-Pitaevskii equation is given by the product u. In 
dimensionless units we define this as the nonlinearity constant Cnl = utolfixi. 



' This restriction simply allows us to avoid using cumbersome notation to account for dilferent spectral bases in each 
direction, and the fully anisotropic case is of no additional computational complexity. 
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In dimensionless units, the PGPE takes the form 

= Ho^c + 'Pc {CnlI<AcI'<Ac) , (A2) 

where 

Ho--\v^ + ^(x^ + y^+f). (A3) 
The c-field is expanded in a spectral basis as 

iAc(x, t) = Yj Cn{t) 0„(x), (A4) 

neC 

where {0„(x)} are the harmonic oscillator eigenstates of with respective eigenval- 
ues e„, and the {c„} are complex amplitudes. The projection is expUcitly implemented 
by limiting the summation indices in ( |A4[ | to the set of values 

C = {n:en< fcut}, (A5) 
i.e. the field i^c only contains the modes of interest. 



A.3 Mode evolution 

Having used the modes of as the spectral basis and to realize the projector, we 



follow the Galerkin approach (i.e. projecting ( A2 1 on to our spectral basis) to obtain 
the amplitude evolution equation 

dcn 

= -/ [e„c„ + CnlG„] , (A6) 

ot 

where 

G„^ J^/^x^:(x)|<Ac(x,?)P<Ac(x,r), (A7) 

is the nonlinear matrix element. Once this matrix elements is evaluated, the evolution 
of the system can be calculated using numerical algorithms for systems of ordinary 
differential equations, e.g. the Runge-Kutta algorithm (e.g. see [171]). Since this is 
a well-understood area of numerical mathematics we do not concern ourselves with 



the details of the propagation algorithm, but instead focus on evaluating ( A7 1. 
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We can point out the central issue for numerical implementation. Expanding the 



fields in expression ( A7 1 into the mode basis we obtain 



Gn = Yj{[ 4>:m;mgimriv] c;c,cr. (as) 

pqr > 

While the matrix elements within the brackets can be exactly calculated in advance, 
computing all G„ values using this expression requires 0{M'^) floating point op- 
erations, where M is the number of c-field region modes. Such scaling would be 
prohibitive for performing realistic calculations. In what follows we show how to 
compute these matrix elements with a scheme that only requires 0{M'^l^) operations. 
Such spectral representations have also been considered for the zero-temperature 
(non-projected) Gross-Pitaevskii equation in references [11,56]. 



A.4 Separability 

An important feature of the basis states (i.e. eigenstates of ^o) is that they are sepa- 
rable into ID eigenstates, i.e. 



0„(X) ifia(x)(fi/j{y)(Py(z), 

e„ <r^ £a + £fi + £y, 
Cyi ^ Co-jSy) 



(A9) 
(AlO) 
(All) 



where {(fa{x)} are eigenstates of the ID harmonic oscillator Hamiltonian, i.e. 



1 cf 1^2 



(A12) 



with eigenvalue e„ = (a -i- ^), for a a non-negative integer. 
For clarity we use Greek subscripts to label the ID eigenstates, so that the specifi- 



cation of the c-field region in ( A5 1 becomes 



C = {a,j3, y:£a + £f} + £y< fcut}- 



(A13) 



Within the c-field region there exists {~ e^ud distinct ID eigenstates (i.e. ipa) in 
each direction, and thus M ~ ^M], 3D basis states (0„) in the c-field region (see the 



left subplot of figure All. 
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A.5 Evaluating the mati-ix elements 

An important observation made in reference [56] was that the nonUnear matrix el- 



ement given in (A7) can be computed exactly with an appropriately chosen Gauss- 
Hermite quadrature. To show this we note that because the harmonic oscillator states 
are of the form ^a(jc) = haHa{x) exp{-x^ /2), where Hdx) is a Hermite polynomial 
of degree a, the field (at any instant of time) can be written as 

.Ac(x, t) = Q{x, y, ~z)e-^''^y'^'''^'\ (A14) 

where 

Q(x,y,z)^ ^ Cafij{T)haHc,(x)hf3H/}(y)hyHy{z), (A15) 

{a/}-y\€C 

is a polynomial that, as a result of the cutoff, is of maximum degree - 1 in the 
independent variables. 

Similarly, it follows that because the interaction term ( A7 1 is fourth order in the 
field, it can be written in the form 

Gafiy = jd'xe-^^'"^'^'^'"^P,f,y(x,y,~z), (A16) 

where 

Pai3y{x,y,z) = haHa(x)h/jHfj{y)hyHy(z)\Q(x,y,z)fQ(x,y,z), (A17) 

is a polynomial of maximum degree 4(M;f - 1) in the independent variables. To 
evaluate these integrals, we note the general form of the Nq point Gauss-Hermite 
quadrature 

X+oo 
dxw(x)m^Yj"'jf{xj), (A18) 
;=i 

where w{x) is a Gaussian weight function, and the Nq values of Wj and Xj are the 
quadrature weights and roots, respectively (see [3]). This quadrature is exact if f(x) 
a polynomial of maximum degree 2Nq - 1 . 



Identifying the exponential term in ( A16 1 as the usual weight function for quadra 



ture, the integral can be exactly evaluated using a three-dimensional spatial grid of 
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8 {Mjf - 1)^ points (i.e. 2 (M^- - 1) points in each direction[^, i.e. 

Ga/3y ^ ^ WiWjWkPa/3y{Xi, Xj, Xk), (A19) 
ijk 

where jc, and are the 2 (Mjc-l) roots and weights of the ID Gauss-Hermite quadra- 
ture with weight function w{x) = exp(-2jc^) [3]. Note, the due to the isotropy of the 
trapping potential, the quadrature grids in all spatial directions are identical. 

A.6 Overview of numerical procedure 

'^{x) quadrature grid representation of field 
c spectral representation of field _■ G nonlinear matrix elements 




Fi gure Al . Schematic of numerical procedure to evaluate the nonlinear matrix elements Gafiy- 



Here we briefly overview how the quadrature described above can be efficiently 
implemented numerically. We require the transformation matrices, given by ID basis 
states evaluated on the quadrature grid, i.e. 

Via - f>a{Xi), (A20) 

to be pre-calculated. Starting from the basis set representation of the field (i.e. {Capy]) 
at an instant of time t, the steps for calculating the matrix elements are as follows 
(also see figure [AT] ): 

(i) Transform from spectral to spatial representation: 

lAc(x,7/t, ?) ^ UiaUjfiUkyCafiy{t), (A21) 



' Since a polynomial of degree 2N - 1 is integrated exactly using an A'-point quadrature. 
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where xijk ^ (xi, xj, xj,). 



(ii) The quadrature integrand of the nonlinear matrix element (A7i is constructed 
by appropriately dividing by the weight function and pre-multiplying by the 
weight^ i.e. 

giXijk) = WiWjWke^^*'''^' \^c{%jk, t)\^^c{%jk, t). (A22) 

(iii) Inverse transforming these integrand functions yields the desired matrix ele- 
ments: 

Gapy = UlU]i,Ul^g{%jk). (A23) 

ijk 

The slowest step in this procedure is carrying out the basis transformation, which 
requires 0{M\), i.e. 0(M^^-^) floating point operations when carried out as a series 
of matrix multiplications. Typical simulations, where we evolve a c-field field with 
M X 2, 000 modes for 100 trap periods, take about 2 hours. 



Appendix B: Numerical technique for the uniform system 
B.l Spectral representation 

The basic quadrature arguments presented for the harmonic oscillator case can be 
applied to the numerical description of the uniform Bose gas. In this section we will 
briefly discuss the uniform case, referring to results from Appendix [A] where they 
are the same. The system of interest is taken to be in a cuboid volume with linear 
dimensions {L, L, L} and subject to periodic boundary conditions. 



The dimensionless PGPE takes the same form as in equation ( A2 1, but with a basis 
Hamiltonian of the form 



{2ny 

where we take periodic boundary conditions, 

>l'cix+ l,y,z) = il^cix,y + 1,2) = ij/cix,y,z + 1) = ^c{x,y,z), (B2) 
and have used xq = L and tQ = mL^/nh as the units of length and time. 



'Here we form e^l^O'l as this corresponds to the polynomial (P) required for the quadrature (see A19 ). 
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As in the harmonic case (see equations (|A9[) - (|A11[)) the basis states are separable 



into ID eigenstates (i.e. ^„(x) = ipa{x)ipi}{y)ipy{z)) of the form 

ip,{x) = e^-\ (B3) 



with the wavevectors ka chosen as harmonics of the perodicity interval, i.e., ka = 
2na, with a an integer, and respective eigenvalues Sa = a^. The values of the indices 



{a,jS, y} specifying the c-tield region is given by equation (A13l, which defines a 
sphere of radius V^cut in a^Sy-space for the uniform system. For later convenience 
we define ffmax as the maximum value of a that occurs in C, i.e. the highest order 
basis state in each direction. For the planewave case we have amax - V?cut> and thus 
in the c-field region we have a total of = 2amax + 1 distinct ID basis states (i.e., 
ifa) in each direction, with M w 3D basis states (i.e., 0„). 



B.2 Evaluating the matrix elements 



In the planewave spectral representation the PGPE takes the form (A6l, for which 



the main challenge is evaluating the nonlinear matrix element (A7). We now show 
how a quadrature approach can be used to evaluate this matrix element exactly. The 
essence of this approach is to transform the field to a spatial representation where the 
nonlinear term is local. 

In each spatial dimension, the quadrature grid of interest (for the uniform case) 
consists of A^Q equally-spaced points given by 

Xj = jAx, l<j<NQ, (B4) 

with spacing Ax = 1/Nq, which spans the spatial region (0, 1]. The quadrature ex- 
pression for an integral of an arbitrary function / is 



pi 

I dxw(x)f{x)xY,wjf{xj), (B5) 
^0 P 



where w(x) = 1 is the weight function, and wj = Ax. That is, for the planewave 
approach, the quadrature rule is the well-known rectangular rule from elementary 
numerical analysis. 

The requirement that our quadrature will exactly calculate the nonlinear matrix el- 
ements is equivalent to the requirement that the ID integrals between are all products 
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of four ipj{x) are evaluated exactly, i.e. 



hpyS = ^ ^X(pl{Xj)(p*p{Xj)ipy{Xj)^ps{Xj), -ffmax < a;,/3,y,6 < ffmax, (B6) 

(B7) 



which holds for the quadrature described above if we take Nq > 2M,,. Thus the most 
efficient and accurate representation is when we choose 2M^ grid points in each 
spatial dimension. 



B.2.1 Fourier interpretation 

The quadrature grid requirement (A^q = 2M^) can be interpreted in terms of Fourier 
properties of the spatial grid. To represent a maximum wavevector of ^cut = 27rQfmax * 
nMx, the Nyquist requirement for the spatial grid is that the distance between points 
should be Ax - l/Mv (or smaller), which requires atleast points over the interval 
(0, 1]. However, our quadrature argument above was that to evaluate the nonlinear 
matrix elements correctly we need twice as many grid points, i.e A^q = 2M^. Such 
a grid is sufficient to satisfy the Nyquist condition for wavevectors of magnitude up 
to 2^cut- To understand why we need so many points consider the worst case for 
the matrix element in ( |B6| ): the case -a = -ji = y - 6 = a^m, i.e. whe re al l 
modes occur with the maximum magnitude wavevector ^cut- The integrand of (B6 1, 
i.e. the product of these four modes, is itself a planewave with wavevector 4^cut. On 
the spatial grid with IM^ points this cannot be represented unambiguously (i.e. it 
exceeds the Nyquist limit of 2^cut)> and is aliased. However, for the choice of 2M^- 
points, this aliasing does not map the wavevector into the region [-^cut> ^cut]> and 
hence does not effect the matrix elements evaluated for the c-field region. For any 
fewer points the aliased wavevector maps into [-^cut. ^cut]> and gives rise to spurious 
dynamics. 



B.3 Overview of numerical procedure 



We could apply an identical procedure to that discussed in section A.6 to evaluate 
the matrix elements with a computational cost per evaluation of 0{Ar). However, 
for the planewave case the basis transformation between spectral (momentum space) 
and position space quadrature grids (i.e. steps (i) and (iii) in section [AT6] ) is equivalent 
to a fast Fourier transformation, which has a computational cost of O [m^ log(M)). 
For more details on the planewave procedure we refer the reader to reference [23]. 
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Appendix C: Mapping to stochastic equations 

The utility of phase space methods requires that the equation of motion for the quasi- 
probabiUty distribution (here ( [58| )) can be mapped to an equivalent stochastic differ- 
ential equation, which is comparatively much easier to solve. A projected functional 
Fokker-Planck equation of the form 



5P 

~dt 



r d\{- ^-^A ((Ac(x), <Ac(x), t) + h.c. 

+ , , Dn («Ac(x), <Ac(x), t) + h.c. 
5(/'c(x)<5iAc(x) ^ ' 

6^ 



-Dn (<Ac(x), (Ac(x), t) + h.c.jp, (CI) 



5(Ac(x)<5iA^(x 



with drift vector A = [A, A*] and diffusion matrix D = [Dn, D12; D*2, D* J has 
an equivalent stochastic equation if the diffusion matrix is positive semi-definite. A 
factorization of the diffusion matrix in the form D = BB^ can then be found, and the 
stochastic differential equation is given by 



rf(^c(x, t) = Pc {A(x, t)dt + B{x, t)dW{x, t)} (C2) 



where di/fc = [difrc,difr*f,]^ , Pc = ['Pc-,'P*qV and dW(x,t) is a vector of noises. 
In general, mapping to ordinary stochastic differential equations is only possible if 
the equation of motion for the quasi-probability is strictly a Fokker-Planck equation 
(derivatives up to second order). 



There is an important technical point regarding the equivalence of (CI I and (C2i. 



The strict equivalence holds only if the projector Pq is implemented with sufficient 



care. In the language of section A.5 the quadrature chosen to compute (C2) must 
be sufficient to generate a c-field delta function of the appropriate numerical order 
upon stochastic averaging. The standard proofs of equivalence [72] can be adapted 
to show that for a diffusion term of polynomial degree 2Dv, giving B of degree D^, 
the delta function must to be a true delta function for terms up to order -1- M^, 
requiring expansion of the noise up to states of degree D^+M,^ and implementation of 



equation ( C2 1 using a numerical quadrature rule sufficient to integrate terms of order 



2Dx + 2Mj^ or a rule of order Z);^ -1- - 1 to generate the appropriate equivalence. 
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